Rare jackpot individuals drive rapid adaptation in Threespine Stickleback

Domain: biology
Published: 2026-03-30
Editorial outcome: major_revision

Abstract
Recombination has long been considered the primary mechanism to bring beneficial alleles together, which can increase the speed of adaptation from standing genetic variation. Recombination is fundamental to the transporter hypothesis proposed to explain precise parallel adaptation in Threespine Stickleback. We study an instance of freshwater adaptation in the Threespine Stickleback system using whole genome data from an evolutionary time-series to observe the genomic dynamics underlying rapid parallel adaptation. Here, we show that rapid adaptation to a freshwater environment depends on a few individuals with large haploblocks of freshwater-adaptive alleles (jackpot carriers) present among the anadromous founders at low frequencies. Biological kinship analyses indicates that mating among jackpot carriers and between jackpot carriers and non-jackpot individuals led to an increase in freshwater-adaptive alleles within the first few generations. This process allowed the population to overcome a substantial bottleneck likely caused by the low fitness of first-generation stickleback possessing a few freshwater-adaptive alleles born in the lake. Additionally, we find evidence that the genetic load that emerged from population growth after the bottleneck may have been reduced through an increase in homozygosity by inbreeding, ultimately purging deleterious alleles. Recombination likely played a limited role in this case of very rapid adaptation.

Introduction
Darwin originally proposed that evolution by natural selection is a gradual process whereby changes in phenotypes that led to speciation occur by selection of variants with small fitness effects over millions of years1. However, there are now numerous examples of rapid evolution, i.e., within tens of generations2–5, which provide an opportunity to study the molecular basis of evolution in more complex eukaryotes outside of a lab setting.

Rapid adaptation proceeds predominantly from standing genetic variation (SGV)6–10, and Gasterosteus aculeatus (i.e., Threespine Stickleback fish) is a premier example of this phenomenon11. In this species, hundreds of loci, most estimated to be more than a million years old12, are the basis for much more recent, precisely repeated, rapid adaptation of ancestral marine or anadromous (i.e., collectively ‘oceanic’) ecotypes to freshwater habitats12–14. This process can occur within a decade12,15.

According to the transporter hypothesis16 (see also ref 10), ancient SGV is maintained at low frequencies in oceanic populations through persistent gene flow from multiple freshwater populations, which are sympatric with anadromous populations during the breeding season. Oceanic-freshwater hybrids with haploblocks (large chromosomal segments with multiple contiguous loci) of freshwater-adaptive alleles move back into the oceanic populations, where they should be disfavored but survive well enough to backcross to anadromous stickleback. Successive generations of backcrossing with oceanic stickleback and recombination break these blocks of freshwater-adaptive alleles apart, and most oceanic individuals are observed to possess only a small number of them, mostly in the heterozygous state12. Consequently, adaptation to new freshwater habitats is thought to proceed from colonizing oceanic individuals, each possessing only a few freshwater-adaptive alleles. In freshwater habitats, selection then favors the freshwater-adaptive alleles, leading to an increase in their frequencies, and genetic linkage among these alleles will allow for their re-assembly back into large haploblocks over multiple generations.

Alternatively, rapid freshwater adaptation may involve selection that strongly favors jackpot carriers. Jackpot carriers are rare oceanic individuals that possess large haploblocks of freshwater-adaptive alleles. Previous studies suggest that the frequency of jackpot carriers in marine environments is about 0.1%12,17. Bassham et al.17 proposed that rapid freshwater adaptation may depend on such jackpot carriers. In addition, simulations have found that a few individuals with large haploblocks of freshwater-adaptive alleles within founding oceanic populations would be sufficient to adapt to new freshwater habitats18. The presence of such individuals could dramatically increase the speed with which haplotypes of adaptive alleles spread through the population compared to the time required for recombination to bring these adaptive alleles together within individuals. In addition, previous work has shown that recombination is significantly suppressed in regions of the genome containing freshwater-adaptive alleles, allowing the maintenance of what has been termed “adaptive cassettes”12,19,20.

In this study, we examine a recently founded freshwater population in Scout Lake, Alaska, derived from anadromous ancestors21. We sequence hundreds of whole genomes from samples collected during the first few generations of rapid adaptation to conditions in Scout Lake. Our results show that adaptation to the freshwater environment depends on the presence of a few jackpot carriers among the founders. After a bottleneck in the third year after founding, subsequent population growth was driven primarily by breeding among the descendants of those few jackpot carriers, leaving signals of marked inbreeding. The increase in homozygosity from inbreeding facilitated the removal of deleterious alleles, thereby likely improving the population's adaptability. Our empirical study shows that only a few jackpot carriers can drive rapid adaptation of anadromous Threespine Stickleback to new freshwater habitats. These findings also suggest a limited role of recombination during the earliest stages of rapid adaptation, in contrast with predictions from the transporter hypothesis.

Results
The Scout Lake experiment
Scout Lake (60.5353N, 150.8322W) is on the Kenai Peninsula, Alaska, USA. After treatment with rotenone in 2009 to exterminate an invasive fish species, 3047 anadromous Threespine Stickleback from a spawning run in Rabbit Slough (62.6893N, 149.427E) in the Matanuska-Susitna Borough, Alaska, were released into the lake in the summer of 2011 (Fig. 1A). Tests with trout late in the summer of 2009 indicated that fish did not survive the rotenone treatment(Supplementary Note 1). The lake is isolated from other bodies of water, and sticklebacks (and other fish) were not detected by exhaustive sampling after the rotenone poisoning. Stickleback samples have been collected annually from this lake since the year after its founding in 2011 (i.e. from 2012) for phenotypic and genetic analysis12,21. In this study, we focused on samples collected in 2013, 2014, 2015, 2017 and 2020 and we referred to these samples as SC2013, SC2014, SC2015, SC2017 and SC2020, respectively. These samples correspond to 2, 3, 4, 6 and 9 years after the founding of the new Scout Lake population. In addition, we included two samples from the anadromous source population, Rabbit Slough, collected in 2009 and 2019. These two timepoints from the Rabbit Slough represent the genomic diversity in the ancestral population over the same period as examined in Scout Lake. We denote these two Rabbit Slough samples as RS2009 and RS2019. Despite being collected ten years apart, the two Rabbit Slough samples showed highly correlated allele frequencies (Fig. S1 and S2) and highly consistent results when used as baselines for inferring processes in the ancestral oceanic population.

The adults stocked in Scout Lake in 2011 produced offspring after release but did not survive the winter (Fig. 1B and Fig. 3 from Kurz et al.22). The F1 progeny likely did not successfully reproduce until 2013. However, small sexually mature males and females were captured in 2012 (Fig. 1B,1C, and 1D; Kurtz et al.22). This finding indicates that the F1 survived in the lake for two years, which is not surprising given that F1 offspring of anadromous Threespine Stickleback typically spend part of the earliest and reproductive segments of their life cycle in freshwater23,24 and thus are physiologically capable of tolerating freshwater. These phenotypic results show that SC2014 represents one-year-old F2 offspring of the F1 generation (Fig. 1E), as any stickleback born in the sampling year would have been too small to be captured with the sampling methods employed (i.e. minnow traps). Thus, the first generation of offspring of the parents that were born and spent their entire life in the lake were adults captured in 2014.

Jackpot carriers increased in frequency in Scout Lake
To infer whether haplotypes with numerous freshwater-adaptive alleles were assembled after introduction to the lake or were already present in rare large haploblocks in a few of the introduced anadromous founders, we used a TN5 transposase-based approach to sequence 432 genomes at coverages ranging from 0.43X to 1.84X (low-coverage). There were 336 low-coverage genomes from individuals collected from two (n = 96), three (n = 48), four (n = 96) and six (n = 96) years after the Scout Lake population was founded and 96 genomes from RS2019. We sequenced additional 20 genomes from SC2020 at a mean of 25.7X coverage (high-coverage) (Methods, Supplementary Data 1), and a lso included 20 previously published high-coverage genomes from RS200912.

We defined multi-SNP haplotypes for each of 344 loci with freshwater-adaptive alleles that we previously identified to experience rapid and significant frequency increases in three lake stickleback populations from Cook Inlet, Alaska12,25. These populations, which included the Scout Lake population, were founded by anadromous ancestors. The 344 loci range from a few base pairs to kilobases, with a median size of 27.3 kb. Each multi-SNP haplotype consists of the most significant freshwater-adaptive SNP within the locus and selected neighboring SNPs whose freshwater allele frequencies are highly correlated (r> =0.99) across our time-series data from the three Threespine Stickleback populations12. We excluded all loci with less than three SNPs, leaving 300 haplotypes that contain three to 3658 SNPs.

We then developed an approximate genotype likelihood approach suitable to call the diploid state of freshwater-adaptive loci as homozygous oceanic, heterozygous, or homozygous freshwater from low-coverage genomic data (Methods). We validated this method by applying it to experimental crosses between parents with known oceanic and freshwater ecotypes derived from a lake (60.914517°, -149.101545°) at mile 87 of Seward Highway a few kilometers east of Girdwood, Alaska (Methods). Genotypic states at loci inferred using our likelihood method were consistent with expectations from our crosses (Fig. S3). Parents were scored as homozygous for the alleles associated with their ecotype at most loci. More importantly, the offspring of parents with alternative homozygous states for a locus were heterozygous at the called locus (Methods, Fig. S3). We then applied this approach to the whole genomes from Rabbit Slough (2009 and 2019) and Scout (2013 to 2020). Loci with missing genotype calls in any individual were excluded, resulting in 280 out of the original 344 loci where all individuals in the study had their genotypes called.

Freshwater-adaptive loci were almost always homozygous for oceanic alleles in the samples from the ancestral population and SC2013 (2 years since founding) (Fig. 1F). Individuals in these samples were either completely homozygous for oceanic alleles (0 out of 20 RS2009; 12 and 15 out of 96 specimens for RS2019 and SC2013, respectively) or possessed freshwater-adaptive alleles at only a few loci (average number of freshwater-adaptive alleles in RS2019 = 2.1 [range: 0-9], SC2013 = 2.25 [range: 0-30], Fig. 1F). The mean proportion of freshwater adaptive alleles in the 20 RS2009 individuals were slightly higher than in either the RS2019 and SC2013 samples (RS2009 = 0.0122 vs RS2019 = 0.00569: t=5.7224, p=0.000008, RS2009 vs SC2013= 0.005952: t=4.9428, p=0.000020). The mean proportion of freshwater-adaptive alleles did not differ significantly between the RS2019 and SC2013 samples (RS2019 = 0.00569 vs SC2013= 0.005952: t=-0.3475, p=0.728697; Fig. 1F, Fig. 2A) or RS2009 and SC2013 samples (RS2009 =0.000875, SC2013= 0.005952: t= –1.7967, p=0.0732, Fig. 2A).

We estimated the size of a contiguous block of freshwater-adaptive alleles by the genome coordinates that defined SNPs. For a series of adjacent loci with freshwater-adaptive alleles, we estimated the size of the block by calculating the difference between the first SNP of the first locus with a freshwater-adaptive allele and the last SNP of that series. For the results presented below, we did not impose a maximum threshold distance between adjacent loci when defining a haploblock. However, we tested the impact of imposing a maximum threshold and this criteria did not change the overall interpretation of the haploblock distributions (Methods, Fig. S4-S5). Of the 81 individuals in SC2013 with freshwater-adaptive alleles, only seven had two or more loci forming a contiguous block with freshwater-adaptive alleles, with blocks of alleles spanning between 0.17Mb to 3.86Mb [0.41cM to 4.57cM] (Fig. 3A, Supplementary Data 2).

In stark contrast to the SC2013 sample, in the third year after founding (SC2014, F2 generation, Fig. 1E), nearly half of the individuals in our sample possessed large numbers of freshwater-adaptive alleles (Fig. 1F, 2B, Fig. 2A). The mean proportion of freshwater-adaptive alleles in SC2014 was significantly different from that of SC2013 (mean proportion of freshwater alleles in SC2014=0.1222 vs SC2013=0.0048; t = 6.3156, p = 8.923e-08; Fig. 3B, Fig. 2A). The distribution of the number of freshwater-adaptive alleles per individual in the SC2014 sample was bimodal (Fig. 2A), which suggested that two distinct groups were in the SC2014 sample (Fig. S6). One group of individuals appear to possess a few number of freshwater-adaptive alleles similar to the ancestral population and the SC2013 samples, while the other had large numbers of freshwater-adaptive alleles (Fig. 2A).

Therefore, to aid our understanding of genomic change during adaptation to freshwater, we categorized individuals as jackpot carriers and non-jackpot individuals using Bayesian Gaussian Mixture Modeling (BGMM)(Methods, Fig. 2B). The BGMM categorized individuals with 5% freshwater adaptive alleles as jackpot carriers with a posterior probability of 0.95 and those with 8% freshwater adaptive alleles were classified as jackpot carriers with 100% certainty (Fig. S7). Based on the BGMM classification, there were no jackpot carriers in RS2009 and RS2019 samples and one jackpot carrier in SC2013 sample (i.e., 1%). There were 23 out of 48 in SC2014 (i.e., 48%), 84 out of 96 in SC2015 (i.e., 88%), and all in SC2017 and SC2020 (i.e., 100%, Supplementary Data 1). We observed that 27 out of the 47 individuals in SC2014 had contiguous blocks with freshwater-adaptive alleles. These contiguous blocks were mostly heterozygous (Fig. S8) and spanned 0.07Mb [0.044cM] to 22.33Mb [70.48cM], with a mean of 2.40Mb [2.9cM] (Fig. 3A and 3B, Supplementary Data 2).

Threespine Stickleback in the Cook Inlet lakes become reproductively mature at around one to two years old,26 and results from Kurz et al.22 indicated that by 2014, the captured individuals were F2 (Fig. 1E). Therefore, the observed jackpot carriers are likely not products of re-assembly via recombination from anadromous founders which possessed only a few freshwater-adaptive alleles. Jackpot carriers must have been present in the founding population, but at such low frequencies that they were undetectable in our limited sampling from the ancestral population (RS2009 and RS2019) or in the SC2013 sample. Based on previous estimates of the frequency of individuals that can be classified as jackpot carriers in ancestral marine environments of 0.1%12,17, the probability of including at least one jackpot carrier in a sample of 96 individuals is only 9% (assuming binomial sampling probabilities). The rapid increase in the frequency of jackpot carriers from 1% in the F1 generation (i.e., SC2013), which was born and spent its entire life in the lake, to almost 48% in the F2 generation (i.e., SC2014) was likely a result of their markedly greater fitness in the freshwater environment compared to individuals with few freshwater-adaptive alleles.

Low survival and/or low reproductive success of the vast majority of non-jackpot carriers during years two and three after founding resulted in a population bottleneck that was manifested by a markedly reduced catch per unit effort (CPUE), a relative measure of abundance (Fig. 4A). In the SC2013 sample, the decline in CPUE was consistent with a decrease in genome-wide genetic diversity compared to the ancestral population samples (Watterson’s Θ=0.005328 for RS2009 and Θ=0.00535 for RS2019 vs 0.0048 for SC2013, Fig. 4B). CPUE declined further in SC2014, and a correspondingly more significant reduction in genetic diversity of almost 50% compared to the ancestral population (Θ = 0.00359 for SC2014, Fig. 4B). This population decline in Scout Lake caused Tajima’s D to become increasingly positive in the SC2013 (D=0.263) and the SC2014 (D=1.618) sample despite the ancestral population having a negative value (D=-0.722 for RS2019) that likely reflects a population expansion in the ancestral environment. The genome-wide site frequency spectrum (SFS), estimated using ANGSD27, indicates progressive, genome-wide loss of singletons, with a particularly drastic distortion of the SFS in SC2014 for rare alleles and is consistent with the dramatic population decline within the Scout Lake population (Fig. 4C). The SFS between the high-coverage RS2009 and SC2020 (Fig. S9) samples also showed loss of singletons in the Scout Lake persisting through several generations after founding.

After this bottleneck, the CPUE increased, and there was a relative increase in the frequency of rarer variants in the SFS for the SC2015 and SC2017 samples compared to SC2014. The number of freshwater-adaptive alleles increased consistently over time in samples collected after the bottleneck; the mean proportion of freshwater-adaptive alleles in SC2015, SC2017 and SC2020, were 0.2078, 0.3504, 0.4292, respectively (Fig. 1F, 2B, Fig. 2A). However, the size of freshwater haploblocks, as measured by the haploblock length in physical base pairs or centiMorgans (Fig. 3), did not change at the same pace, as would have been expected if recombination had been bringing freshwater-adaptive alleles together from multiple parents each generation to form increasingly larger haploblocks.

There was no significant difference in the mean haploblock lengths between the two Rabbit Slough samples (RS2009 vs RS2019; t(158.42) = 1.2021, p = 2.311e-01 ), as well as between the RS2019 and SC2013 (t(334.40) = 0.4335, p = 6.645 e-01; Fig. 3, S4, S5). There was a statistically significant difference between the mean haploblock lengths of SC2013 and SC2014 (t(802.17) = -9.7188, p = 3.4958e-21), SC2014 and SC2015 samples (t(1515.87) = -2.2363, p = 2.5478e-02) and SC2015 and SC2017(t(4475.41) = -2.7869, p = 5.3432e-03). There was no statistically significant difference between the mean haploblock lengths of SC2017 and SC2020 (t(748.07) = -1.4509, p = 1.4724e-01, Fig. 3). We further explored the differences in haploblock lengths using Cohen’s D, which expresses the mean difference relative to the pooled standard deviation. We estimated Cohen’s D and bootstrapped confidence interval for each pair of consecutive timepoints. There was a small difference between RS2009 and RS2019 (D=-0.1585[-0.352 – 0.087] ) and RS2019 and SC2013 (D= 0.0044[-0.1341 – 0.26] ) compared to SC2013 and SC2014 (D=0.4025[0.367 – 0.4510] ). Between SC2014 and SC2015, there was a small difference in the mean haploblock lengths (D=0.093[0.012 – 0.172]), suggesting that despite the statistically significant t-test, the difference was not substantial at this time period. We observed similar Cohen’s D between mean haploblock means between SC2015 and SC2017 (D=0.083[0.024 – 0.138]) as well as SC2017 and SC2020 (D=0.074[-0.025 – 0.174]). These results suggest that there was limited increase in the haploblock lengths across the time periods, especially after the apparent demographic bottleneck in 2014, which would be expected if the haploblocks were being reassembled by recombination.

Venu et al.19 indicated that as many as half of all chromosomes in male, as well as a third of chromosomes in females, could be inherited without any recombination event per generation. In addition, recombination is suppressed in individuals that are heterozygous for marine and freshwater haplotypes at adaptive loci. Our sampling period covered nine years, spanning 5-6 generations, and the majority of the individuals sampled after the bottleneck were predominantly heterozygous at freshwater adaptive loci. We estimated linkage disequilibrium (LD) across freshwater adaptive loci (Methods) to assess whether recombination had any substantial effect on the integrity of haploblocks. We observed an elevated LD from SC2014 decay compared to SC2013 and samples from the ancestral population, consistent with a strong directional positive selection increasing the frequency of adaptive haploblocks or potentially the effect of demographic bottleneck that occurred in 2014 (Fig. 5). In subsequent years, we found no observable patterns of LD decay, suggesting that there was not enough time for recombination to have substantial effect on the integrity of haploblocks.

The homozygosity of freshwater-adaptive loci also increased over time from SC2014 (Fig. S8), which could reflect the mating between heterozygous jackpot carriers. These results demonstrate that haploblocks of freshwater-adaptive alleles were not re-assembled through progressive recombination over many generations, but instead the frequencies of freshwater-adaptive alleles increased after the F1 generation by matings predominantly between the descendants of jackpot carriers present in the founders. These jackpot carriers likely had an increased fitness in the freshwater environment.

Biological kinship and inbreeding during freshwater adaptation
Our sequence data from multiple whole genomes during the earlier years after the Scout Lake population was founded provided a unique opportunity to explore the genealogical relationships among individuals as the population adapted to the new freshwater environment and to infer the dynamics of haploblock inheritance. We estimated biological relatedness among our low-coverage genomes using READv228,29 and ngsRelate30. The estimates of relatedness coefficients from the two methods were highly correlated (r>0.92 for all time points; Fig. S10). However, in situations in which there is increased allelic sharing compared to expectations from panmixia due to inbreeding, we might expect READv2 relatedness coefficients to be more robust as the approach relies on observed rates of allele mismatching, while ngsRelate utilizes allele frequencies under the assumption of Hardy-Weinberg equilibrium, leading to potentially inflated values of relatedness (Methods). Therefore, we report the results based on READv2, which identifies up to third-degree relatives. All pairwise relationships detected using READv2 and NGSrelate are provided in Supplementary Data 4.

We also estimated individual inbreeding coefficients, F, with ngsRelate, which can estimate inbreeding from low-coverage genomes (up to approximately 4X30). A down-sampling experiment of the high-coverage genomes from the SC2020 sample demonstrated that ngsRelate underestimates inbreeding coefficients from low-coverage of 1X by a factor of ~50% compared to high-coverage (~30X) whole genome data (Methods, Fig. S11); thus the reported inbreeding coefficients from our low-coverage genomes from the SC2013 to SC2017 samples are likely underestimated.

There were no close biological kinships(up to third-degree relatedness) in samples taken from the ancestral population, suggesting that the founders of the Scout Lake population were drawn from an outbred anadromous population with a large effective size (Ne). Despite possessing identical proportions of freshwater-adaptive alleles as the ancestral population (Fig. 1F), the SC2013 sample had two first-, 38 second- and five third-degree relatives (Fig. 6A, S12A), indicating that this population is descended from a small number of fish with a much more limited choice of mates than the ancestral anadromous population.

Some related pairs involved one individual from the SC2013 sample and the other from the SC2014 sample. There were 36 third-degree and two second-degree related pairs between these two time points (Fig. 6A). Some related pairs had one of the pair being a jackpot carrier in SC2014 and the other being a non-jackpot individual in SC2013, suggesting that jackpot carriers mated with non-jackpot individuals early on in the colonization process.

While the majority of relatedness estimates within the SC2013 sample (Fig. S12A) formed small networks, each containing only a few individuals, out of a total of 60 related pairs within 2014, 57 were connected within a single large network (Fig. S12B). The one first-degree relationship was a sibling pair. Apart from four relationships between jackpot and non-jackpot individuals, this network's remaining pairs were all jackpot carriers (Fig. 6A, S12B). Most jackpot carriers (20 out of 23 individuals) identified in the SC2014 sample were connected within this single network. The third-degree relationships of jackpot carriers sampled in 2014 are likely to be cousins or half-avuncular assuming the typical lifespan of freshwater Alaskan Threespine Stickleback is one to two years31.

We observed one first-degree (sibling pair), nine second-degree, and 315 third-degree pairs of relatives between specimens from the SC2014 and SC2015 samples. Within the SC2015 sample (Fig. S12C), there was one first-degree (sibling pair), 36 second- and 298 third-degree pairs. These 335 connections were between 77 out of 96 sampled individuals in SC2015, 55 of which were jackpot carriers (there were 78 jackpot carriers in SC2015). Most individuals sampled in SC2014 and SC2015 were members of the large network of relatedness which contained most of the sampled jackpot carriers (20 out of 23 in SC2014 and 67 out of 78 jackpot carriers in SC2015; Fig. 6B). Unlike in the SC2013 and SC2014 samples, some individuals in SC2015 showed evidence of inbreeding (Fig. 6D), with F values ranging from 0.0 to 0.149 (mean=0.0289). Thus, the jackpot carriers observed in SC2015 were mostly related to the SC2014 jackpot carriers, suggesting that the SC2015 jackpot carriers descended from the network of jackpot carriers in SC2014.

The network of relationships between the SC2015 and SC2017 sample provided further evidence that population growth was driven by jackpot carriers, with 34 second- and 1201 third-degree pairs in 167 out of 192 (i.e. 87%) individuals from SC2015 and SC2017. All relationships, except seven third-degree and one second-degree relative, were between jackpot carriers. By 2017, all sampled individuals were jackpot carriers, and 93 out of them were in a single network (Fig. S12D). The individuals sampled in SC2017 had an average inbreeding coefficient of 0.0316 (range [0.0, 0.1121]). None of the 20 individuals in the SC2020 sample were first-, second- or third-degree relatives, though we identified three third-degree relatives between SC2017 and SC2020 individuals (Supplementary Data 4). To test if this lack of relatedness within SC2020 was significantly lower than in SC2017 or simply a result of the smaller sample size (SC2017 n = 96, vs. SC2020 n=20), we performed a permutation test by randomly sampling 20 individuals from the SC2017 sampled 1,000 times and observing the number of pairwise relatives of at least third-degree. In no permutation did we observe 20 random individuals from SC2017 with no relatedness (mean=57.2, min=16, max=112, Fig. S13). This result suggests that the lack of relatedness in the SC2020 sample was likely to be a result of population growth that reduced the probability of sampling relatives and not a result of the small sample size. Despite the deficiency of related pairs, individuals in SC2020 were, on average, inbred with an F of 0.078 (range [0.038, 0.110]). The dearth of relatives within SC2020 suggests that the population had become large enough to allow more mate choice by this time, with the increased F reflecting the legacy of inbreeding from the preceding generations.

Patterns of polymorphism, selection and genetic load
When a population experiences a decline in size, the increased probability of inbreeding can increase genetic load32–34 and reduce mean population fitness34,35. Rare, deleterious, recessive alleles previously masked in heterozygotes become exposed as homozygotes and reduce individual fitness32,33,36. To determine how the dynamics of the observed bottleneck affected the genetic load of the Scout Lake population, we estimated the site frequency spectrum (SFS) genome-wide as well as at 4-fold (silent sites) and 0-fold (amino acid replacement) degenerate sites for samples from the time series (Methods).

The proportion of rare alleles in an SFS provides insight into the demographic history of a population, with excess rare alleles suggesting population growth, while a deficit may suggest population decline. This pattern results from each individual introducing unique mutations to the population as the population grows. Demographic processes affect both 0- and 4-fold sites equally, whereas natural selection is more likely to impact 0-fold sites. As expected, 0-fold sites generally demonstrate approximately a third of the diversity of 4-fold and genome-wide sites for all samples, presumably due to the effects of purifying selection (Fig. 7H). The ancestral population (RS2019) showed a slight excess of singletons at 0-fold sites compared to 4-fold and genome-wide sites, possibly due to the segregation of slightly deleterious recessive alleles in the ancestral anadromous population at low frequencies (Fig. 7A, S14A). During the bottleneck from 2013 to 2014 (SC2013 and SC2014 samples), there was a similar reduction in Watterson’s Θ at 0-fold and 4-fold sites and genome-wide (30%, 34% and 25%, respectively; Fig. 7D-E), indicating that a demographic factor, the population bottleneck, shaped most of the patterns of diversity during this time (Fig. 7H).

However, within the SC2014 sample a notable increase in singletons at 0-fold sites relative to that expected at genome-wide sites was detected, with 4-fold sites being intermediate, presumably reflecting linkage to 0-fold sites (Fig. 7D). Compared to non-jackpot individuals, the SFS of jackpot carriers was skewed towards more common variants with no significant increase in 0-fold singletons, likely reflecting the close genealogical relationships between these individuals (Fig. S15). The singletons in non-jackpot carriers reflected the excess of singletons at 0-fold sites in the adult F2 generation (SC2014). We could, therefore, attribute the observed singletons at 0-fold sites almost entirely to non-jackpot carriers. However, this excess of singletons at 0-fold sites was not observed in the samples collected from 2015, 2017 and 2020 (SC2015, SC2017 and SC2020). As described above, we also observed high levels of relatedness and inbreeding in SC2015 and SC2017 as shown by increased homozygosity (Fig. S8). Thus, slightly deleterious 0-fold sites were likely being purged from the population as a result of a) non-jackpot individuals failing to leave descendants except when they mated with jackpot carriers and b) an increased chance of recessive deleterious alleles being homozygous and thus being selected against.

Testing jackpot-mediated rapid adaptation with forward-in-time simulations
To further explore whether our observed data at Scout Lake could be better explained by the inheritance of large haploblocks of freshwater adaptive alleles, rather than their gradual accumulation through recombination, we employed forward-in-time Wright-Fisher simulations in SLiM37 to model the process of freshwater adaptation. Our simulations closely follow the scenarios in Galloway et al. 202018 and Roberts Kingman et al. 202112 (Supplementary Note 4, Fig. S16), with the latter using a machine learning inference framework to identify a selective regime that accurately fits observed allele frequency trajectories at adaptive loci across multiple rapidly adapting freshwater populations, including Scout. Based on the demographic and selective parameters inferred in Roberts Kingman et al. 2021, we let the simulation run for 1000 generations, which was enough time for the system to reach equilibrium, and then founded a new freshwater population under two scenarios: A) a scenario where both jackpot individuals (maintaining the definition of 5% freshwater alleles) and non-jackpot carriers are present in the founding population of Scout Lake and B) a scenario where only non-jackpot carriers are present in the founding population of Scout Lake.

We estimated the total freshwater adaptive alleles in the founding population under both scenarios (Fig. 8A), and found that the freshwater content was not significantly different (KS test: 0.08, p-value:0.908). Thus the presence of a small number of jackpot individuals in the founding population (Fig. S17) did not skew the overall content of freshwater alleles available for adaptation for scenario A compared to B. We then calculated the proportion of freshwater adaptive alleles in the generation corresponding to the observed sampled timepoints from Scout Lake. We found that we could only replicate the observed mean proportions of freshwater adaptive alleles in simulations under scenario A (Fig. 8B), demonstrating that jackpot individuals are needed in the founding population for the observed rapid adaptation in Scout Lake. When jackpot carriers were removed from the founding population, rapid adaptation was not possible in any simulation (Fig. 8C). Thus our simulations indicate that it is crucial to have individuals with large blocks of freshwater adaptive alleles for rapid freshwater adaptation as we observed in Scout Lake.

Methods
Sample collection, library preparation and genome sequencing: All fish used in this study were collected according to an approved protocol from the Institutional Animal Care and Use Committee (IACUC 1446584) at Stony Brook University and from The College of New Jersey Institutional Care and Use of Animals Committee (protocols 1908-001MW1A1 and 2002-001MW1A3). Samples of Threespine Stickleback have been collected from Scout Lake at least once a year around the start of the breeding season in late May or June since 2012. Minnow traps with a mesh of 6.35 mm or usually 3.175 mm have been set yearly for up to 24 hours at less than 2 m depth and 5 m from shore. After capture, sticklebacks are separated from other fish species and transferred to a bucket of lake water to which equal volumes of sodium bicarbonate and MS-222 (tricaine methane-sulfonate) had been added at a high enough concentration to cause the fish to lose equilibrium within 30 seconds and to die within a few minutes. Death was inferred from the failure of the fish to react to tapping the side of the bucket in which they were held and then to pinching the caudal fins of selected fish. We netted the fish out of the bucket, washed them in lake water, and dropped them into 70 % ethanol in deionized water in a 1-liter bottle with up to about half fish and the remainder with ethanol solution. The ethanol was replaced with a fresh 70% ethanol solution within about 24 hours after the lipid from the fish had discolored the ethanol with a yellow hue. The right pectoral fin was usually clipped from the specimen and placed into a small, numbered conical tube of 70% ethanol for DNA extraction, and the remainder of the fish was placed with another fish from the same sample in a numbered, 15-ml tube with 70% ethanol and one fish head up and the other head down. These fish are stored either in David M. Kingsley’s laboratory in the Department of Developmental Biology, Stanford School of Medicine or in the Veeramah lab at Stony Brook University.

DNA was extracted using the Qiagen DNeasy 96 Blood & Tissue Kit for animal tissue and was quantified with a Qubit Fluorometer 3.0 set to High Sensitivity option for dsDNA. For the high-coverage DNA sequences from SC2020, we selected high quality samples with minimal evidence of fragmentation to sequence. The average concentration of the 20 selected samples was 72.32 ng/uL, ranging from 41.2 ng/uL to 134.8 ng/uL. We sent the DNA for resequencing at the Beijing Genomics Institute (BGI) using their proprietary DNBseq technology. For the low-coverage genomes, we randomly selected DNA from 96 specimens with high concentrations of DNA for each time point (except SC2014, in which there were only 48 specimens) for library preparation with plexWell™ 384 (SeqWell, Beverly, MA, USA). The plexWell™ 384 kit is based on seqWell’s proprietary TN5 transposase that involves the insertion of Illumina i7 adapters into individual DNA samples (up to 96 samples), which are then pooled together. An Illumina i5 adapter is inserted into each pool and amplified in a PCR reaction using manufacturer-provided library primer mix and Roche’s KAPA HiFi DNA Polymerase readymix. The concentration of each pool was within the manufacturer-recommended range of 4-8ng/μL (RS2019-4.2 ng/μL, SC2013-5.68ng/μL, SC2014-6.84ng/μL, SC2015-4.56ng/μL, SC2017-5.24ng/μL, Crosses-7.56ng/μL). We performed further quality controls by running aliquots of the prepared libraries on 2% agarose gel and an Agilent 2100 Bioanalyzer. We pooled all individuals on a 96-well plate together and added a unique i5-barcoded adapter. Each population of pooled 96 individuals is then sequenced on an Illumina HiseqXten platform, with paired ends read length of 150 bps at an output of around 45Gb. Threespine Stickleback has a genome size 460Mb, so 45Gb output yields approximately 100X coverage.

Crosses from Mile 87 lake for validation of called genotypes at freshwater adaptive loci: In early June of 2021, adult sticklebacks were collected from Mile 87 Lake (60.914517°, -149.101545°) using unbaited, 10mm mesh minnow traps set within a few meters of the shoreline. Anadromous and freshwater males and females in breeding conditions were transported live in aerated lake water to the temporary laboratory facility in nearby Anchorage, Alaska. We identified ecotype identity based on substantial size and body shape differences between the adults and the anadromous (larger) and freshwater (smaller). We produced four crosses via in vitro fertilization: male anadromous x female anadromous, male freshwater x female freshwater, male anadromous x female freshwater, and male freshwater x female anadromous. We euthanized males in an overdose of buffered MS-222, and their testes were dissected and macerated in a drop of embryo medium (2ppt Instant Ocean Sea Salt in distilled water + 1 drop of methylene blue/2L). We stripped ovulated eggs from a gravid female into a Petri dish and subsequently covered them with the extracted sperm. After three minutes to allow fertilization, we rinsed clutches thoroughly, submerged them in fresh embryo medium, and stored them at 7oC. The following day, we separated embryos, removed dead/unfertilized individuals, and refreshed the embryo medium. We maintained clutches each day. Prior to five days post fertilization, clutches were shipped overnight to The College of New Jersey in Ewing, NJ. Embryos began hatching ~12-14 days pf, and fry was reared for an additional three weeks when they were euthanized in an overdose of buffered MS-222 and preserved in 100% ethanol for subsequent DNA extraction. Fish were collected under Alaska Department of Fish and Game permit P-21-006.

Sexual maturity scoring of Scout 2012 and 2013 samples: We removed the gonads of females and males to perform macroscopic inspections of them. Ovaries of females were assigned to six reproductive phases following57,58: latent, early maturing, late maturing, mature, ripening, and ripe. These phase assignments have been updated with more recent terminology based upon microscopic histological observations59 but were used as initially designated. Latent females were considered sexually immature, whereas females in the other five phases were considered sexually mature. Testes of males were assigned to latent or mature phases based on their size and opacity. Latent males had small, transparent or translucent testes, and mature males had enlarged, cloudy to opaque white testes.

Catch per unit effort estimation: Catch per unit effort (CPUE) is the total number of fish caught divided by “trap hours” The number of trap hours is the product of the number of traps and the number of hours we set. We captured all Threespine Sticklebacks using 1/8 or ¼-inch mesh Gee Minnow Traps set overnight for <24 hours within <5 m of shore and in water <1.5 m deep. We placed the traps in and adjacent to patches of rooted aquatic plants if they had begun to grow. The bottom was mostly sand or silt. We set about 100 traps behind the same three residences on the north shore of the lake.

Bioinformatics Processing. Sequencing reads were de-multiplexed into individual genomes based on the i7-barcodes by BGI. The total number of bases sequenced per sample ranged from 0.004 Gb to 1.79 Gb, with a mean of 0.54 Gb. The total number of reads per sample ranged from 21,266 to 11,917,120, with a mean of 3,028,575 reads. After trimming adapters with AdapterRemoval (ver. 2.2.2)60, we mapped reads to the Threespine Stickleback genome version gasAcu1-412,61. Libraries sequenced across multiple Illumina lanes and runs were processed separately, and we added read groups using Picard before being merged with samtools 62. We also marked duplicates using Picardtools. We performed base recalibration using BaseRecalibrator from GATK version 3.763,64.

Genotype calling and validation. As mean coverage was ~1x across samples, we developed a novel method for calling diploid genotype states at freshwater-adaptive loci. As described in the Results, we identified highly linked SNPs within loci that underlie freshwater adaptation, using the correlation coefficient (>0.99) of SNP allele frequencies from Pool-Seq experiments covering a time series of lake populations undergoing freshwater adaptation: Loberg, Cheney and Scout Lakes12 (Supplementary Data 5). These multi-SNP haplotypes allowed us to call the diploid state of each locus based on an approximate genotype likelihood approach that employed multi-SNP haplotype information. Genotype likelihoods were calculated for individual SNPs within each locus using the approach described in DePristo et al.65 and then summed across correlated SNPs based on the freshwater or anadromous allele to obtain approximate likelihoods for anadromous homozygous, heterozygous and freshwater homozygous allelic states for each locus. Genotype calls were assigned for a locus if there was a minimum of 3 SNPs with reads with mapping error probabilities less than 0.001 and base call error probabilities less than 0.01. We coded freshwater-adaptive loci that did not meet these filtering criteria as missing.

We validated the ability of this method to determine the diploid state of the loci by applying it to call haplotypes of a set of experimental crosses between anadromous and freshwater parents from the impoundment at mile 87 along Seward Highway. We observed that estimated states from our likelihood method were highly consistent with expectations from the results of our crosses (Fig. S1). Parents were scored as homozygous for the alleles associated at most loci with their ecotype. More importantly, offspring in these experimental crosses inherited states consistent with the states of their two parents. We did not consult that parental genotypes were to call genotypes of their progeny. For example, the offspring of parents with alternative homozygous states at the locus were invariably scored as heterozygous.

We also imputed our low-coverage data with beagle 4.066 to enhance genotyping accuracy using a reference panel of 169 medium- to high-coverage whole genomes (see below). We then adjusted our likelihood approach to incorporate imputation-based genotype probabilities rather than individual SNP genotype-likelihoods. The overall pattern of the genotypic state of each locus was not qualitatively different between our likelihood approach and Beagle-imputed genotype probabilities (Fig. S1), as the dosage of freshwater-adaptive alleles carried by each individual was similar in both methods. Beagle imputation led to a reduced rate of missingness per locus (280 loci without missingness after vs. 235 before imputation with Beagle).

Imputation of low-coverage data with Beagle 4.0: Using beagle 4.0 and a reference panel of 169 genomes, we imputed the missing data in the low-coverage sequences. The reference panel comprised genomes mapped to the Gas-Acu1-4 genome, with bwa-mem and base recalibration performed on the bams with GATK’s BaseRecalibrator using hard-filtered SNPs from Roberts-Kingman et al.12. The coverage of the included genomes ranged from 9–62X (Supplementary Data 3). It included 69 oceanic fish, 58 from well-established freshwater populations, 40 from experimental transplants, and two of unknown ecotype designation with a global distribution from 12,67. We followed GATK’s best practices guidelines for variant calling with HaplotypeCaller, joint genotyping, and variant quality score recalibration with VariantRecalibrator. We only included biallelic sites in the reference panel and targeted 6,027,376 high-quality SNPs. We downloaded beagle 4.0 from https://faculty.washington.edu/browning/beagle/beagle.r1399.jar to perform haplotype phasing per chromosome and imputation on all low-coverage data simultaneously. The genotype likelihood used was estimated using GATK’s Haplotypecaller. We used recombination maps from ancestral Rabbit Slough, which was previously estimated12 using LDhelmet68. We filtered the recombination maps to exclude absent sites in the reference panel. Then we ran beagle 4.0 as follows java -Xmx90G -jar beagle.27Jan18.7e1.jar gl=chrII.gz ref=WGS_170_recalibrated_Snps6mil_chrII.phased.vcf.gz map=chrII.map impute=True out=All_samples_chrII. After imputation, we combined all chromosomes using bcftools concat.

Haploblock construction and erosion: To assess the impacts of imposing a gap threshold on the observed distribution of defined haploblocks for the freshwater adaptive loci (Fig. 3, supplementary Data 2), we analysed the haploblock distributions while imposing a physical distance threshold. We first filtered sites with missing data (as was the case for our initial analyses presented in the main manuscript). We chose the maximum gap threshold as twice the size of the largest adaptive locus, which was ~0.4Mb (excluding the three well-known inversions on chrI, XI and XXI, which are significantly larger than most adaptive loci and span up to 1.9Mb). Thus, if adjacent loci with freshwater alleles are more than 0.8Mb apart, we break any single larger haploblock into two smaller haploblocks. When we compared the haploblock distributions when imposing a threshold of 0.8Mb (Fig. S4, S5) with those estimated in the original analysis, we do see smaller haploblocks in both physical and genetic distances, which is to be expected as blocks are broken up. However, more importantly, when imposing this 0.8Mb threshold we still see no major increase in haploblock size across time points using both physical and genetic distance compared to the rapid increase in per individual freshwater content, consistent with our original inference from Fig 3, and there is a very strong correlation in mean haploblock size with and without imposing the 0.8Mb threshold (Fig. S18).

Classification of individuals as jackpot carriers: We performed a dip test to check the modality of the SC2014 sample based on freshwater content (i.e. number of freshwater alleles per individual) and identified two modes. We also pooled all samples across timepoints for the same metric of the number of freshwater alleles per individual, and used a Bayesian Gaussian Mixture Modelling (BGMM) approach implemented in scikit-learn (version 1.7.1) in python 3.11.13 to estimate the number of underlying distributions that best fit the data. BGMM incorporates a Dirichlet process prior that allows the model to infer the effective number of components directly from the data. We set the maximum number of components to four and used the full covariance type, which treats each component to have its own general covariance matrix. However, doubling the number of maximum components did not change the result of the analysis. Bayesian information criterion (BIC) suggested two components to be the best choice, consistent with the dip test (Fig. S6). The fitted model provided posterior probabilities for each observed individual’s membership in each of the two main components based on their freshwater content. We then categorized all the datapoints at each timepoint into one of these two components (Fig. 2). The model categorized individuals with 5% freshwater adaptive alleles as jackpot carriers with a posterior probability of 0.95. Individuals with 8% freshwater adaptive alleles were classified as jackpot carriers with 100% certainty (Fig. S7).

Estimating the site frequency spectrum (SFS) with theta and Tajima’s D. To account for our low-coverage genome data, we used RealSFS from ANGSD69 to estimate each sample's site frequency spectrum (SFS) based on genotype likelihoods. For SC2015, we identified four duplicates through our relatedness analyses, which were removed before estimating the SFS. We masked ChrXIX and all transposons before estimating the SFS. We used the ancestral genome constructed by Roberts Kingman et al.12. Then we imported the estimated SFS into dadi70 and folded the SFS, and performed further analyses, including summary statistics like nucleotide diversity (π), Watterson’s estimator of Θ and Tajima’s D71. To prevent any bias introduced by unequal sample sizes of various time points, we used dadi to project down each time point to 40 samples. For 0-fold and 4-fold degenerate sites, we used degenotate (https://github.com/harvardinformatics/degenotate) to compute the degeneracy of coding sites across the stickleback genome. In total, we identified 24,891,530 0-fold and 6,348,969 4-fold degenerate sites. Then, we used the -sites option from ANGSD to estimate the SFS for the 0-fold and 4-fold sites.

Linkage disequilibrium (LD) decay: To estimate the LD decay across the sampling period, we focused on the 280 adaptive loci without any missing data. For each time point, we estimated pairwise r2 among all combinations of the 280 loci using hard-called locus-based genotypes, coded as 0 for homozygous freshwater, 1 for heterozygous, and 2 for homozygous marine, corresponding to the number of marine alleles carried by each individual at each locus. To quantify the uncertainty in LD estimates, we performed bootstrapping using genotype matrices for each of the timepoints. For each timepoint, the genotype matrix was resampled with replacement 1,000 times to generate a distribution of r2 values in order to estimate 95% confidence intervals. We then binned pairs of loci/SNPs by distance to examine differences in the rate of LD decay.

Kinship analysis. We utilized two methods to infer pairwise relatedness, both designed for low-coverage data, READv229 and NGSrelate30. We called genotype likelihoods genome-wide for all SNPs identified previously by Roberts Kingman et al.12 using GATK Unified Genotyper v3.7.065. We estimated allele frequencies for each Scout Lake sample separately using the approach of Kim et al.72. We then intersected these allele frequencies with allele frequencies estimated from pool-seq data from Roberts Kingman et al.12 from the same time points using the approach of Lynch et al.73. For each time point, any SNP that had a total depth of coverage (across all samples) of less than half or double the mean coverage for both the TN5 and pool-seq data was included in the genome mask. In addition, we performed Fisher’s exact test on a 2x2 contingency table containing the depth of each allele in the TN5 and Pool-Seq data and added any SNPs with a p-value less than 0.05 to the mask. We merged masks from each time point, known repeats, and transposable elements. Then, we applied this mask to all time points, resulting in a set of 2,119,671 high-quality SNPs for relatedness analysis. We made haplodiploid genotypes for use in READv2 based on the highest allele depth28,29, which estimates the degree of relatedness between pairs of individuals by using the proportion of allele mismatches across the genome. We used READv229, which estimates up to third-degree relatedness. We validated READv2 results using ngsRelate30, which utilizes population allele frequencies to estimate identity-by-descent probability. We ran ngsRelate with options ngsRelate -F 1 -h input.vcf -O output -z samp_list -A AF, where AF is the allele frequency estimated from a set of unrelated samples from SC2013 (again estimated using the approach of Kim et al.72 )

We used Cotterman coefficients74 for first-degree relatives to distinguish parent-offspring from sibling pairs. We relied on the Cotterman coefficients estimated from ngsRelate instead of READv2. When there is no inbreeding, the Jacquard coefficients J9, J8 and J7 map to the Cotterman coefficients. Based on the Cotterman coefficients, the five first-degree relatives observed in our Scout Lake samples were likely siblings. However, READ2 designates 4 out of these five as parent and offspring. In READv2, the Cotterman coefficients k0 and k2 correspond to the proportions of windows classified as unrelated and identical. This proportion is expected to be low for parent-offspring pairs and higher for siblings.

Inbreeding Estimation. We used NgsRelate30 to estimate individual inbreeding coefficients using allele frequencies from unrelated specimens from SC2013. To validate the estimates of inbreeding coefficients from ngsRelate (Fig. 6D), which is designed for low-coverage datasets, we used PLINK to estimate the runs of homozygosity for our high-coverage SC2020 samples. We compared the results with estimates of inbreeding from ngsRelate for them. We estimated runs of homozygosity by first performing LD pruning in PLINK with options --indep-pairwise 50 5 0.5. We then estimated the runs of homozygosity with PLINK, option --homozyg --homozyg-window-kb 5000. There is a high correlation between the runs of homozygosity and estimated inbreeding coefficients from ngsRelate (r=0.86, Fig. S19).

We also downsampled the high-coverage SC2020 samples to 1X coverage using samtools (v1.9) with option samtools view -s x -b input.bam -o downsampled.bam, where x is calculated as 1/y, and y is the coverage of the high-coverage sample estimated using GATK’s DepthOfCoverage. We used ngsRelate to estimate the inbreeding coefficient using allele frequencies from SC2013. There is a high correlation between estimates from downsampled and high-coverage samples (r=0.822, Fig. S11).

Forward Wright-Fisher simulation: We performed forward-in-time Wright Fisher simulations to model the process of adaptation with and without jackpot individuals using SLiM v537,75. We assumed non-overlapping generations, with a generation time of one year, following the framework of 12 and 18. We simulated a large oceanic (anadromous) population connected to 10 distinct freshwater populations (Fig. S16). We allow migration to occur between the freshwater populations and the marine population, but not among freshwater populations. As in Roberts Kingman et al.(2021), we set the migration rate from marine to freshwater populations (M_AN_TO_FW) to be 0.001, and migration rate from freshwater populations to marine population (M_FW_TO_AN) to be 0.01. At the start of the simulation, we introduce freshwater-adaptive mutations at all 341 loci in one of the freshwater populations to mimic an established population that is fully adapted to a freshwater environment. We do not allow de novo mutations during the simulation and each mutation has a selection coefficient of 0.01 in the freshwater environment. The value of 0.01 is of similar magnitude to the mean selection coefficient inferred across loci in Roberts Kingman et al.(2021). At generation 1000, we established a new freshwater population representing Scout Lake by randomly sampling individuals from the anadromous population. Following its founding, no further migration was permitted between this simulated Scout Lake and any other populations in the system.

Phenotypic characterization. We measured three traits that differ consistently between anadromous and freshwater stickleback76 and have diverged rapidly from the ancestral anadromous condition in other recently found lake populations in the SC2014 sample. The three traits were lateral plate morphs, gill rakers and standard length (Supplementary Note 2).

Statistics and reproducibility:We performed statistical analyses in Python. Sample sizes are reported in the figure legends. All data and code necessary for reproducing the conclusions in this study have been made publicly available.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript by Kwakye et al. addresses an intriguing question about the genomic mechanisms underlying rapid adaptation to freshwater environments in Threespine Stickleback, a model organism for evolutionary biology. This work challenges the traditional recombination-centered \"transporter hypothesis\" by demonstrating empirically that rapid freshwater adaptation can be mediated primarily through rare individuals carrying large haploblocks of adaptive alleles (jackpot carriers), rather than gradual recombination events. The study leverages an exceptionally rich temporal genomic dataset and offers insights into kinship dynamics, population bottlenecks, genetic load management, and the initial absence of phenotypic divergence, enhancing our understanding of evolutionary mechanisms. Nevertheless, several points should be addressed to strengthen the manuscript further:\n1. The authors used a 2019 Rabbit Slough sample as the ancestral reference for detecting standing genetic variation underlying rapid adaptation. However, the actual introduction event occurred in 2011, making a 2009 or 2011 sample more appropriate as the baseline. The choice of a 2019 baseline may introduce biases from genetic drift or selection occurring over the intervening period. The authors should clearly justify using the 2019 sample or preferably reanalyze key results with the available 2009 data to ensure robust identification of adaptive loci.\n2. The authors mention that the SC2013 samples may include a mixture of F0 founders and early-born F1 individuals. However, the manuscript does not explain why no genetic substructure or differentiation between these two groups is observed within the SC2013 cohort. This point should be explicitly addressed.\n3. The authors define jackpot carriers as individuals with more than 10% freshwater-adaptive alleles, based on the bimodal distribution observed in the SC2014 sample. While this threshold is reasonable, it would strengthen the methodological rigor if the authors could briefly discuss whether alternative thresholds (e.g., 8% or 12%) were tested and how sensitive their key results (e.g., haploblock expansion patterns, relatedness networks) are to the exact cutoff chosen. Additionally, it would be informative to comment on whether the definition of jackpot carriers remains consistent across later timepoints (e.g., SC2015, SC2017) as selection proceeds and freshwater dosage distributions shift.\n4. While the manuscript defines a contiguous haploblock as consisting of two or more consecutive freshwater-adaptive loci, it remains unclear whether a maximum physical distance threshold was imposed between adjacent loci, or how missing data across loci were handled during block construction.\n5. Given the progressive increase in inbreeding coefficients across timepoints, it would be valuable for the authors to briefly discuss whether increased homozygosity could have influenced the observed patterns of haploblock expansion or adaptive allele frequency changes.\n6. As recombination could theoretically break apart large freshwater haploblocks over successive generations, it would be useful if the authors could briefly discuss whether any evidence of haploblock erosion (e.g., internal recombination events, partial block loss) was observed by SC2017 or SC2020.\n7. While the authors clearly state that Scout Lake was cleared of predators and competitors through rotenone treatment before stickleback introduction, they do not discuss how this simplified ecological context might have influenced the rapid adaptation process observed. The absence of ecological interactions such as predation and interspecific competition could fundamentally alter selective pressures, potentially making Scout Lake a special case. I recommend that the authors explicitly address this point in the Discussion and temper the generalization of their findings accordingly.", "reviewer_2": "This study by Kwakye and colleagues investigates the genomic mechanisms of rapid freshwater adaptation in Threespine Stickleback. Using time-series genomic data, they propose that adaptation was primarily driven by rare \"jackpot\" individuals present in the founders, who carried large, pre-existing haploblocks of freshwater-adaptive alleles. These individuals and their descendants are suggested to have rapidly increased in frequency following a population bottleneck, leading to inbreeding and potential purging of deleterious alleles. This proposed mechanism offers a valuable contrast to the traditional transporter hypothesis.\nMajor Suggestions:\n1- The central claim regarding the indispensability of jackpot carriers for the observed rapid sweep would be substantially strengthened by quantitative modelling. Employing explicit ABC or forward simulations could test this hypothesis against alternatives and help disentangle selection from genetic drift, particularly during the bottleneck.\n2- the argument for early physiological adaptation preceding morphological changes, currently based largely on negative evidence, would be more compelling with direct positive evidence. Measuring relevant physiological traits, such as osmoregulation or metabolic rates, in early-generation individuals could provide this. Or tune down it or move to the supplemenary to save text in the main text (see below).\n3- the manuscript's length and density could be reduced for clarity and impact. Now it is too long for the journal. Streamlining the narrative to focus on the core \"jackpot carrier\" hypothesis and condensing or moving less central results, such as the detailed morphology and fads2 findings, to supplementary information is recommended.\nMinor Suggestions:\n1- To improve accessibility for a broader audience, the authors should enhance the in-text explanation of Fig. 6, more explicitly guiding readers through the interpretation of SFS patterns and their connection to conclusions about genetic load and selection.\n2- all custom scripts and bioinformatic pipelines used for the genomic analyses should be made available.\n3- The figure legend for Fig. 1F: the 344 freshwater-adaptive loci, should briefly explained.\n4- The Y-axis labels for the two columns in Fig. 2 need to be clarified.\n5- Fig 4 panel labels should be annotated"}, "author_response": "RESPONSE TO REVIEWER COMMENTS FOR:\nRARE JACKPOT INDIVIDUALS DRIVE RAPID ADAPTATION IN THREESPINE\nSTICKLEBACK” (NCOMMS-25-20702-T)\nOur responses to the reviewer’s comments are below and any text added to the main text are\nmarked red.\nReviewer #1 (Remarks to the Author):\nThe manuscript by Kwakye et al. addresses an intriguing question about the genomic\nmechanisms underlying rapid adaptation to freshwater environments in Threespine\nStickleback, a model organism for evolutionary biology. This work challenges the\ntraditional recombination-centered \"transporter hypothesis\" by demonstrating empirically\nthat rapid freshwater adaptation can be mediated primarily through rare individuals\ncarrying large haploblocks of adaptive alleles (jackpot carriers), rather than gradual\nrecombination events. The study leverages an exceptionally rich temporal genomic dataset\nand offers insights into kinship dynamics, population bottlenecks, genetic load\nmanagement, and the initial absence of phenotypic divergence, enhancing our\nunderstanding of evolutionary mechanisms. Nevertheless, several points should\nbe addressed to strengthen the manuscript further:\n1. The authors used a 2019 Rabbit Slough sample as the ancestral reference for detecting\nstanding genetic variation underlying rapid adaptation. However, the actual introduction\nevent occurred in 2011, making a 2009 or 2011 sample more appropriate as the baseline.\nThe choice of a 2019 baseline may introduce biases from genetic drift or selection occurring\nover the intervening period. The authors should clearly justify using the 2019 sample or\npreferably reanalyze key results with the available 2009 data to ensure robust identification\nof adaptive loci.\nWe tested the correlation between the minor allele frequencies estimated from two Rabbit\nSlough samples (a sample of 100 individuals from 2009 sequenced using pool-seq approach in\nRoberts Kingman et. al1 and the 96 low coverage genomes from 2019 used in this study). We\nfound a very high significant correlation between the two timepoints when considering all SNPs\n(Review figure 1) and when we only look at SNPs within regions of the genome with freshwater\nadaptive alleles1 (Review figure 2). This demonstrates that drift has not had a strong effect in the\nintervening years, which is expected given that the Rabbit Slough population is directly\nconnected to the sea and hence has a large effective population size. We also note that the 2019\nRabbit Slough sample was not used to identify adaptive loci. All adaptive loci we used were\nidentified by Kingman et al. 1, which utilized the RS2009 sample along with pooled samples\nfrom Loberg, Cheney and Scout, including the Scout time points reanalyzed using whole\ngenomes in this study. The RS2019 sample was simply used to demonstrate the typical\ndistribution of freshwater alleles in a large anadromous population sample (figure 1F), and a\nsimilar distribution is observed in the smaller sample of 2009 whole genomes and pools from\nKingman et al (Review figure 3)\nWe have now added the sentence below to the manuscript to explicitly state that the two samples\nare very similar [lines 339-342], and added Review figures 1 and 2 below to the supplementary\nfigures (fig. S3 and S4).\n“The allele frequencies at freshwater adaptive loci, estimated from a sample collected in 2009\nbefore the founding of the lake and the RS2019 sample, were highly correlated (fig. S3, S4),\nconsistent with a large and stable anadromous population.”\nReview figure 1: Correlation between minor allele frequencies\nestimated from pool-seq of 100 individuals sampled in Rabbit\nSlough in 2009 and 96 individuals sampled from Rabbit Slough in\n2019 used in this study, which was used in the main manuscript\nusing all SNPs.\nReview figure 2: Correlation between minor allele\nfrequencies estimated from pool-seq of 100 individuals\nsampled in Rabbit Slough in 2009 and 96 individuals\nsampled from Rabbit Slough in 2019 used in this study,\nwhich was used in the main manuscript using SNPs\nwithin freshwater adaptive regions.\nReview figure 3: Genotypes of 20 whole genomes sampled from Rabbit\nSlough in 2009 showing predominantly marine genotypes.\n2. The authors mention that the SC2013 samples may include a mixture of F0 founders and\nearly-born F1 individuals. However, the manuscript does not explain why no genetic\nsubstructure or differentiation between these two groups is observed within the SC2013\ncohort. This point should be explicitly addressed.\nWe apologize if our original language was not clear about this point, but actually, we did not find\nevidence that supported a mixture of F0 founders and early born F1 individuals. We indicated\nthis on lines 154-155, where we state that “The adults stocked in Scout Lake in 2011 produced\noffspring after release but did not survive the winter (fig. 1B and Fig. 3 from Kurz et al.2).” Thus,\nwe found through our own analyses, and that from Kurz et al2 that the population in SC2013\nwere all F1 progeny, not a mixture of F0 founders and early-born F1 individuals (see fig. 1E). In\naddition, our Fig. 1 length information specifically provided context on the reproduction and\ncohorts in the first 4 years, where there was no evidence for a mixture of the founders and F1 in\nthe year after founding and that the F1s only reproduced in the second year and not as one year\nolds. We do note that we used ‘mostly’ on line 220 in the original manuscript in the sentence\nbeginning “Out of the 81 individuals in SC2013 (mostly F1 generation in the lake)...”. Perhaps\nthis is where some of the confusion arose, where we were attempting to denote that, although our\nevidence does not support their presence in SC2013, it is not impossible that there might be F0\nfounders. Therefore we have now removed “(mostly F1 generation in the lake)” from that\nsentence. Perhaps there is another area of the manuscript we accidently gave the impression that\n2013 contained a mix of F0 and F1s that we missed, and if the reviewer could point it out, we\nwould be happy to amend the text appropriately.\n3 [part A]. The authors define jackpot carriers as individuals with more than 10%\nfreshwater-adaptive alleles, based on the bimodal distribution observed in the SC2014\nsample. While this threshold is reasonable, it would strengthen the methodological rigor if\nthe authors could briefly discuss whether alternative thresholds (e.g., 8% or 12%) were\ntested and how sensitive their key results (e.g., haploblock expansion patterns, relatedness\nnetworks) are to the exact cutoff chosen.\nAs suggested, we adjusted the cut off to 8%, and 12% for calling jackpot carriers and reassessed\nour results. Changing the threshold affected the categorization of just one individual in SC2014,\nwith a 12% cut off leading to one less jackpot carrier. In SC2015, an 8% cut-off led to the\nclassification of three non-jackpot individuals as jackpot carriers, while increasing it to 12%\nexcluded 4 jackpot carriers (Review figure 4). There is also one jackpot carrier in SC2017 that\nhad fewer than 12% freshwater adaptive alleles (Review figure 4). Regardless, the pattern of\ngrowth of proportion of jackpot individuals is qualitatively the same, demonstrating that\nalternative thresholds have minimal effects on which specimen gets classified as a jackpot\ncarrier, and therefore their relative distribution across the timepoints. In fact, the definition of\nwhich samples get classified as a jackpot carrier does not influence the framework of any of our\nanalyses from a statistical perspective, except for our site frequency spectrum (SFS) analysis\nbetween jackpot carriers and non-jackpot individuals in the SC2014 sample in fig. S11, which is\nthe only analysis where the two groups are directly compared.\nReview figure 4: Distribution of jackpot carriers over time using different thresholds\n(8, 10, and 12%) as a definition of jackpot carrier\nOf the few individuals that get re-classified based on shifting the cut off, the majority had one or\ntwo connections in our relatedness networks, except for one individual in SC2015. This\nindividual was classified as a jackpot using a 10% cut off but did not meet the 12% cut off as a\njackpot. The individual had multiple relatives, including a first-degree relative, and has been\nhighlighted in Review figure 5. The fact that this individual is part of the larger cluster of related\njackpot individuals indicates that our 10% threshold was probably the most appropriate in\ndefining jackpots, though in general the structure of the relatedness networks is almost identical\nregardless of the threshold used.\nOverall, shifting the threshold up or down does not lead to conclusions largely different from\nwhat we observed from the 10% we originally set in the manuscript. We have added lines 243-\n245 to make this point and also Review figure 5 to the supplementary figures [fig. S5] to point\nout the effect of altering the threshold.\nReviewer figure 5: Relatedness of samples from SC2014, SC2015 and\nSC2017. The highlighted individual, SC-2015-PB-X153-E11 had multiple\nrelatives. \n3 [part B]. Additionally, it would be informative to comment on whether the definition of\njackpot carriers remains consistent across later timepoints (e.g., SC2015, SC2017) as\nselection proceeds and freshwater dosage distributions shift.\nWe note that our definition of jackpot carrier never changes throughout the study. It is always\ndefined as individuals with >=10% freshwater allelic content. At no point does this definition\nchange. We simply ask whether the number of individuals who meet this criteria change across\nthe various time points.\n4. While the manuscript defines a contiguous haploblock as consisting of two or more\nconsecutive freshwater-adaptive loci, it remains unclear whether a maximum physical\ndistance threshold was imposed between adjacent loci, or how missing data across loci were\nhandled during block construction.\nWe thank the reviewer for this important methodological observation. In the initial manuscript\nwe did not impose a maximum threshold between adjacent loci, and we removed all loci with\nmissing data.\nWe have added the following lines 163-166 to the manuscript to state this explicitly:\n“We then applied this approach to the whole genomes from Rabbit Slough and Scout. Loci with\nmissing genotype calls in any individual were excluded, resulting in 235 loci out of the original\n344 that were analysed in subsequent analysis.”\nAnd lines 183-186:\n“For the results presented below, we did not impose a maximum threshold distance between\nadjacent loci when defining a block. However, imposing a maximum threshold did not change\nthe overall interpretation of the haploblock distribution”.\nHowever, we have also considered the impact of this approach and re-constructed the\nhaploblocks by setting different maximum thresholds between adjacent loci (i.e. a gap\nthreshold). We first filtered sites with missing data (as was the case for our initial analyses\npresented in the main manuscript). We chose the maximum gap threshold as twice the size of the\nlargest adaptive locus, which was ~0.4Mb (excluding the three well-known inversions on chrI,\nXI and XXI, which are significantly larger than most adaptive loci and span up to 1.9Mb). Thus,\nif adjacent loci with freshwater alleles are more than 0.8Mb apart, we break any single larger\nhaploblock into two smaller haploblocks. When we compared the haploblock distributions when\nimposing a threshold of 0.8Mb (Review figure 6 and 7) with those estimated in the original\nanalysis, we do see smaller haploblocks in both physical and genetic distances, which is to be\nexpected as blocks are broken up. However, more importantly, when imposing this 0.8Mb\nthreshold we still see no major increase in haploblock size across time points using both physical\nand genetic distance compared to the rapid increase in per individual freshwater content,\nconsistent with our original inference from Fig 2, and there is a very strong correlation in mean\nhaploblock size with and without imposing the 0.8Mb threshold (Review figure 8).\nReview figure 6: Comparison of physical distance for each\ntimepoint for different cut off scenarios. A) When there is\nno threshold (which is reported in the main manuscript) B)\nWhen the threshold is set to 0.8Mb\nReview figure 7: Comparison of genetic distance for each\ntimepoint for different cut off scenarios. A) When there is no\nthreshold (which is reported in the main manuscript) B) When the\nthreshold is set to 0.8Mb\nReview figure 8: Plot of average genetic distance when setting a threshold of 0.80Mb (y-axis)\nand when no threshold is set (x-axis).\n5. Given the progressive increase in inbreeding coefficients across timepoints, it would be\nvaluable for the authors to briefly discuss whether increased homozygosity could have\ninfluenced the observed patterns of haploblock expansion or adaptive allele frequency\nchanges.\nOn lines [441-442] (which was 662-664 of the original text), we suggested that homozygosity\nresulting from inbreeding due to the reduced population size could potentially have led to\npurging of maladapted alleles. This would have increased the fitness of the population emerging\nfrom the bottleneck and thus perhaps provided a base for faster positive selection from the\nfreshwater alleles at adaptive loci (i.e. quicker allele frequency changes). Beyond this, we are not\naware of any other arguments that could be made regarding the link between homozygosity and\nhaploblock expansion/allele frequency increases.\n6. As recombination could theoretically break apart large freshwater haploblocks over\nsuccessive generations, it would be useful if the authors could briefly discuss whether any\nevidence of haploblock erosion (e.g., internal recombination events, partial block loss) was\nobserved by SC2017 or SC2020.\nRecombination may indeed (almost certainly) have broken some freshwater haploblocks apart.\nHowever this signal would have been overwhelmed by the rapid increase in average per-\nindividual freshwater content resulting from increased rate of mating between individuals with\nexisting freshwater haploblocks, which was still continuing in SC2020 as seen by the increased\nhomozygous freshwater alleles at this time point. As far as we are aware, the only way we would\npossibly be able to observe internal recombination events or block loss in our data is if we had\nparent-offspring trios to directly these processes. Unfortunately, despite a lot of pairwise\nrelatedness in our sample, we do not have any trios to fully address the reviewer’s observation.\n7. While the authors clearly state that Scout Lake was cleared of predators and competitors\nthrough rotenone treatment before stickleback introduction, they do not discuss how this\nsimplified ecological context might have influenced the rapid adaptation process observed.\nThe absence of ecological interactions such as predation and interspecific competition\ncould fundamentally alter selective pressures, potentially making Scout Lake a special case.\nI recommend that the authors explicitly address this point in the Discussion and temper the\ngeneralization of their findings accordingly.\nWe completely agree with the Reviewer that Scout Lake may not capture the full ecological\ncontext of a typical case of freshwater adaptation, given that it was cleared of certain competitors\nprior to the founding. It is something we [the authors] discussed a lot during the preparation of\nthe manuscript, and how to best present our specific results. As a consequence throughout the\noriginal manuscript, we attempted to guard against over-generalization of our findings. In the\nAbstract, for instance, we used the phrase “an instance of freshwater adaptation…” [line 65] to\ndenote that our results pertain to Scout Lake. At the beginning of the Discussion, we stated, “We\nstudied the rapid adaptation of an anadromous Threespine Stickleback population to a new\nfreshwater environment…”. And in the conclusion we noted, “Studies of other genomic time\nseries from other Threespine Stickleback will determine how general our results are…” Thus we\nbelieve we have been careful with generalizing our findings, despite realizing the importance of\nthese findings to our understanding of rapid adaptation in threespine stickleback, as well as its\nimplications to evolutionary biology in general.\nHowever, we have now also added the following sentence [lines 532-535] to state this limitation\nof our study more explicitly.\n“While our experimental set up mimicked natural colonizations as much as possible, we note\nthat Scout Lake represented a somewhat simplified ecological context as it was originally\ncleared of certain competitors through rotenone treatment prior to the founding of the\nThreespine Stickleback population, even though some of these were re-stocked. Future work\nutilizing alternative study designs will be essential to determine how generalizable our results\nare likely to be.”\nReviewer #2 (Remarks to the Author):\nThis study by Kwakye and colleagues investigates the genomic mechanisms of rapid\nfreshwater adaptation in the Threespine Stickleback. Using time-series genomic data, they\npropose that adaptation was primarily driven by rare \"jackpot\" individuals present in the\nfounders, who carried large, pre-existing haploblocks of freshwater-adaptive alleles. These\nindividuals and their descendants are suggested to have rapidly increased in frequency\nfollowing a population bottleneck, leading to inbreeding and potential purging of\ndeleterious alleles. This proposed mechanism offers a valuable contrast to the traditional\ntransporter hypothesis.\nMajor Suggestions:\n1- The central claim regarding the indispensability of jackpot carriers for the observed\nrapid sweep would be substantially strengthened by quantitative modelling. Employing\nexplicit ABC or forward simulations could test this hypothesis against alternatives and help\ndisentangle selection from genetic drift, particularly during the bottleneck.\nWe have now included a forward simulation testing two scenarios with and without jackpot\nfounders, which we believe have significantly strengthened our manuscript, and thank the\nreviewer for their suggestion. We performed Wright-Fisher simulations in SLiM mimicking the\ntransporter hypothesis, as in Galloway et al. 20204 and Roberts Kingman et al. 20211. We let the\nset up run for 1000 generations, after which we founded a new freshwater environment,\nrepresenting the Scout Lake population. In one scenario, we found the new lake with both\njackpot carriers and non-jackpot individuals, and in another removed the jackpot carriers. Our\nresults showed that the rapid adaptation that we observed empirically in Scout Lake could only\nbe replicated in the presence of jackpot carriers. Even one jackpot carrier was enough for rapid\nadaptation. We have added a new section in the main paper to describe these simulations and the\nresults [lines 447-473] and also added the simulation setups and details to the supplementary\ninformation [supplementary section 15]\n2- the argument for early physiological adaptation preceding morphological changes,\ncurrently based largely on negative evidence, would be more compelling with direct\npositive evidence. Measuring relevant physiological traits, such as osmoregulation or\nmetabolic rates, in early-generation individuals could provide this. Or tune down it or\nmove to the supplementary to save text in the main text (see below).\nWe agree that the results of this section are not that compelling (though it provides some\nimportant motivation for future studies beyond the scope of this paper). We have now moved the\nsection on morphological phenotypes and the corresponding discussion to the supplementary\ntexts accordingly, which also helps address the comment below regarding article length.\n3- the manuscript's length and density could be reduced for clarity and impact. Now it is\ntoo long for the journal. Streamlining the narrative to focus on the core \"jackpot carrier\"\nhypothesis and condensing or moving less central results, such as the detailed morphology\nand fads2 findings, to supplementary information is recommended.\nAs recommended, we have reduced the length and density of the manuscript. The original\nsubmission had 8005 words (abstract plus main text–introduction, results, discussions and\nconclusion). We have removed and/or moved certain sections to the supplementary. We have\nremoved the sections on “Rapid freshwater adaptation was not sex-specific” and moved\n“Phenotypes of jackpot carriers and non-jackpot individuals in SC2014” to the supplementary.\nWe have also removed some of our discussions on “The consequences of bottlenecks on genetic\nload during rapid adaptation”.\nIn addition to these major changes, we have also moved/removed certain paragraphs as follows:\nMoved to supplementary:\n“An alternative scenario for the presence of recent hybrids between anadromous and freshwater\nstickleback in the Scout Lake samples is that the jackpot carriers were produced by matings\nbetween anadromous founders and very rare freshwater individuals from the original Scout Lake\npopulation that had survived the rotenone treatment of the lake. We tested this hypothesis by\ncreating hybrids from established freshwater populations from the Pacific Coast and Rabbit\nSlough individuals (see Supplementary Section 13). The distribution of freshwater-adaptive\nalleles in these hybrids was 50 to 63% compared to the 12.5% to 49.5% content we found in the\nindividuals we sampled in the SC2014 sample, suggesting that, the jackpot carriers observed in\nSC2014 sample were probably recent oceanic-freshwater hybrids, or their progeny. These were\nlikely to have been present among the Rabbit Slough founders and not hybrids between survivors\nfrom the original freshwater Scout Lake population and released anadromous individuals. This\nlatter scenario is further supported by previous studies that have found jackpot carriers in\nsamples collected from marine environments.”\nRemoved:\n“First-degree relatives can either be parents-offspring pairs or siblings. Second-degree relatives\ninclude grandparent-grandchild pairs and avuncular pairs. Third-degree relatives include first\ncousins, great aunts/uncles, grandniece/nephews, and great-grandparents. We can distinguish\nthe two types of first-degree relatives but not the various types of second or third-degree\nrelatives”\n“The individual with notable contiguous blocks of freshwater-adaptive alleles in the SC2013\nsample was a third-degree relative of an individual in the SC2014 sample with a similar\nproportion of freshwater alleles. No relatives from this pair were found in the SC2015 sample or\nlater, and this lineage presumably died out in 2014 or failed to proliferate enough to be\nrepresented in subsequent samples, possibly because they possessed too few freshwater-adaptive\nalleles. Thus, their freshwater content of ~7% was insufficient to contribute to the growing\npopulation from 2014 and beyond, justifying our criteria of defining a jackpot carrier as\nindividuals with > 10% freshwater-adaptive alleles.”\n“Our kinship analysis demonstrated that jackpot carriers mated with non-jackpot individuals\nearly on in the adaptive process, prior to the bottleneck in 2014. In order to examine the extent\nto which founding non-jackpot individuals contributed freshwater-adaptive alleles to the later\npost-bottleneck population, we identified loci where all 21 jackpot carriers found in the SC2014\nsample were homozygous for the oceanic allele. There were four such unlinked loci. These loci\nincluded an 11Kb locus on chromosome IV, a 2Kb locus on chrIX, and 11Kb and 19Kb loci on\nchrXIX. We then determined if any individual in the subsequent timepoints (SC2015, SC2017 and\nSC2020) had inherited freshwater alleles at these loci. We found only five heterozygous\nindividuals across all four loci in all the subsequent samples (N=212). However, there were no\nfreshwater adaptive alleles in the non-jackpot individuals observed in the SC2013 and SC2014\nsamples at these four loci. This indicates that non-jackpot individuals within the founders\ncontributed minimally to the pool of adaptive alleles circulating after the bottleneck.”\nThese adjustments together with other editing have resulted in reduction of the manuscript length\nto 5977 words in the main text (which now meets the word limit for the journal) compared to the\noriginal 8005.\nMinor Suggestions:\n1- To improve accessibility for a broader audience, the authors should enhance the in-text\nexplanation of Fig. 6, more explicitly guiding readers through the interpretation of SFS\npatterns and their connection to conclusions about genetic load and selection.\nWe have added the following sentences to lines [413-417] to guide readers on the SFS\ninterpretations we make:\n“In a SFS, the proportion of rare alleles can reflect the demographic history of a population,\nwith excess rare alleles suggesting population growth, while a deficit may suggest population\ndecline. This pattern results from each individual introducing unique mutations to the population\nas it grows, and vice versa during decline. Demographic processes are also likely to affect both\n0- and 4-fold sites equally, while natural selection is more likely to impact 0-fold sites.”\n2- all custom scripts and bioinformatic pipelines used for the genomic analyses should be\nmade available.\nAll custom scripts have been deposited to the first author’s GitHub and can be accessed through\nthe following link: https://github.com/a-kwakye/Rare-Jackpot-Individuals-Drive-Rapid-\nAdaptation-in-Threespine-Stickleback. We have also deposited these scripts and intermediate\ndatasets that can be used to replicate the figures as a capsule on Code Ocean and can be accessed\nonce the manuscript is published https://codeocean.com/capsule/2139373/tree\n3- The figure legend for Fig. 1F: the 344 freshwater-adaptive loci, should briefly explained.\nWe have added a short description of the 344 freshwater-adaptive loci to the legend of figure 1.\n4- The Y-axis labels for the two columns in Fig. 2 need to be clarified.\nWe have updated the Y-axis label to the sample names instead of year since founding to be\nconsistent with the text and other figures.\n5- Fig 4 panel labels should be annotated\nThis has been corrected.\nReferences\n1. Roberts Kingman, G. A. et al. Predicting future from past: The genomic basis of recurrent\nand rapid stickleback evolution. Sci Adv 7, (2021).\n2. Kurz, M. L., Heins, D. C., Bell, M. A. & von Hippel, F. A. Shifts in life-history traits of two\nintroduced populations of threespine stickleback. Evol. Ecol. Res. 17, 225–242 (2016).\n3. Hendry, A. P. et al. Designing eco-evolutionary experiments for restoration projects:\nOpportunities and constraints revealed during stickleback introductions. Ecol. Evol. 14,\ne11503 (2024).\n4. Galloway, J., Cresko, W. A. & Ralph, P. A Few Stickleback Suffice for the Transport of\nAlleles to New Lakes. G3 10, 505–514 (2020)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have made a commendable effort, including adding new analyses (e.g., sensitivity testing of haploblock thresholds and alternate definitions of “jackpot carriers”), and expanding the discussion in several key areas. These additions improve the clarity and scope of the manuscript. However, several methodological concerns remain unresolved, and important interpretive and theoretical issues continue to undermine the strength of the main conclusions. In particular, the current analytical framework still relies on arbitrarily defined thresholds without sufficient theoretical justification or statistical validation. Moreover, potential confounding factors, such as inbreeding-driven homozygosity and recombination-driven haploblock erosion, are not thoroughly examined or incorporated into the interpretation of adaptive dynamics.\nMajor concerns:\n1. The authors argue that the RS2019 sample was not used to identify adaptive loci and that allele frequencies are highly correlated between RS2009 and RS2019. However, this response does not fully address the original concern. The key issue is not the origin of adaptive loci per se, but whether the RS2019 sample can reliably serve as a reference for standing variation prior to the 2011 introduction event. The authors demonstrate that adaptive alleles underwent selection between 2011 and 2015, raising the possibility that some allele frequency changes also occurred in the RS population over the same period. Given this context, using RS2019 to represent ancestral allele frequencies without any temporal correction or bootstrapping introduces potential biases, even if correlations appear high. More importantly, the RS2009 sample is explicitly available (as acknowledged), and its allele frequencies were used in previous publications. There is no compelling reason not to replicate key frequency-based visualizations (e.g., Fig. 1F) using RS2009 data, or at least provide results based on both datasets to confirm robustness. Without this, the choice of RS2019 remains methodologically suboptimal and weakens the interpretation of standing genetic variation.\n2. The classification of “jackpot carriers” using a fixed 10% threshold for freshwater-adaptive alleles is central to multiple downstream analyses, yet the rationale for this cutoff remains insufficiently justified. While the authors briefly explored sensitivity to alternative thresholds (8% and 12%), these appear arbitrarily selected and do not constitute a systematic evaluation. Moreover, Figure 2C shows that the freshwater content distribution shifts substantially across years, suggesting that a fixed threshold may capture biologically distinct subsets of individuals at different timepoints. This raises concerns about the comparability and interpretation of “jackpot” status across temporal scales. The authors are encouraged to provide a stronger conceptual and statistical justification for the 10% threshold, and to consider whether a relative or quantile-based definition (e.g., top x%) might yield more robust or interpretable results given the temporal dynamics of freshwater allele accumulation.\n3. The authors acknowledge that elevated homozygosity resulting from inbreeding may have facilitated the purging of maladaptive alleles, however, this response does not fully address the potential mechanistic implications of increased homozygosity on haploblock structure and allele frequency dynamics. For instance, greater homozygosity could reduce effective recombination, leading to the apparent expansion or preservation of extended haplotypes over time. Additionally, in small populations, hitchhiking effects may be more pronounced, allowing non-adaptive alleles to rise in frequency due to linkage with adaptive loci. These possibilities warrant at least brief acknowledgment in the discussion to avoid over-simplification of the observed genomic patterns as purely selection-driven. A more nuanced consideration of demographic and genetic drift effects would strengthen the interpretation of haploblock dynamics over time.\n4. While the authors reasonably argue that increasing homozygosity could obscure recombination signals, the potential erosion of large freshwater haploblocks remains an important consideration. Even in the absence of parent-offspring trios, population-level patterns, such as reduced haploblock lengths, internal heterogeneity, or linkage disequilibrium decay, could provide indirect but informative evidence of recombination. Currently, the manuscript presents no empirical analysis to support the conclusion that recombination played a limited role. Given that recombination is a central process influencing haplotype integrity and adaptive allele assembly, a lack of data-driven support weakens the strength of this claim. A brief analysis or at least an expanded discussion of potential haploblock erosion would strengthen the interpretation that haploblock persistence reflects biological processes (e.g., selection, demographic dynamics) rather than methodological limitations or overlooked recombination events.", "reviewer_2": "I am satisfied with the authors’ thorough and thoughtful revisions. The incorporation of forward simulations, restructuring of the manuscript for clarity, and relocation of sections to the supplementary material have significantly strengthened the paper. The authors have addressed all of my concerns.\nCongratulations to the authors on an excellent study."}, "author_response": "RESPONSE TO ADDITIONAL REVIEWER COMMENTS FOR:\nRARE JACKPOT INDIVIDUALS DRIVE RAPID ADAPTATION IN THREESPINE\nSTICKLEBACK” (NCOMMS-25-20702B)\nWe want to thank the reviewer for their additional comments, which we believe have\ntremendously improved our manuscript. Our responses to the reviewer’s comments can be found\nbelow and any text added to the manuscript are marked red.\n1. The authors argue that the RS2019 sample was not used to identify adaptive loci and\nthat allele frequencies are highly correlated between RS2009 and RS2019. However, this\nresponse does not fully address the original concern. The key issue is not the origin of\nadaptive loci per se, but whether the RS2019 sample can reliably serve as a reference for\nstanding variation prior to the 2011 introduction event. The authors demonstrate that\nadaptive alleles underwent selection between 2011 and 2015, raising the possibility that\nsome allele frequency changes also occurred in the RS population over the same period.\nGiven this context, using RS2019 to represent ancestral allele frequencies without any\ntemporal correction or bootstrapping introduces potential biases, even if correlations\nappear high. More importantly, the RS2009 sample is explicitly available (as\nacknowledged), and its allele frequencies were used in previous publications. There is no\ncompelling reason not to replicate key frequency-based visualizations (e.g., Fig. 1F) using\nRS2009 data, or at least provide results based on both datasets to confirm robustness.\nWithout this, the choice of RS2019 remains methodologically suboptimal and weakens the\ninterpretation of standing genetic variation.\nAs suggested by the reviewer, we have included the RS2009 sample to Fig. 1F (Review Figure\n1), Fig. 2 (Review Figure 2) and Fig.5 (Review Figure 3). We have also added the RS2009\nsample to the estimated site frequency spectra (Review Figure 4) as supplementary figure S7.\nWe additionally estimated the relatedness between the individuals in the RS2009 sample and the\nlevel of inbreeding. As observed in RS2019, we found no related individuals and no evidence of\ninbreeding in the RS2009 sample. The number of freshwater adaptive alleles in RS2009 ranged\nfrom 0.4% to 2% while RS2019 ranged 0% to 1.3% (Review Figure 5). In addition, all the\nindividuals from both RS2009 and RS2019 are assigned to the non-jackpot category by our new\nBGMM (see below). Overall, both RS2009 and the RS2019 samples appear to be a good\nreference for standing variation in the ancestral population, with the RS2019 having an\nadvantage of having a much larger sample size for downstream analysis such as assessing the\ndegree of relatedness and quantifying the very low frequency of jackpot individuals present in\nthe marine population.\nWe now include RS2009 in all our data interpretation regarding the ancestral marine population\nand relevant figures and have also added the paragraph below at the beginning of the results to\nindicate we use two Rabbit Slough samples spanning the sampling period in Scout Lake.\nLines 110- 113 “In addition, we included two samples from the anadromous source population,\nRabbit Slough, collected in 2009 and 2019. These two timepoints from the Rabbit Slough\nrepresent the genomic diversity in the ancestral population over the same period as examined in\nScout Lake. We denote these two Rabbit Slough samples as RS2009 and RS2019. Despite being\ncollected ten years apart, the two Rabbit Slough samples showed highly correlated allele\nfrequencies (Fig. S2 and S3) and highly consistent results when used as baselines for inferring\nprocesses in the ancestral oceanic population.”\nReview Figure 1\nReview Figure 2\nReview Figure 3\nReview Figure 4\nReview Figure 5\n2. The classification of “jackpot carriers” using a fixed 10% threshold for freshwater-\nadaptive alleles is central to multiple downstream analyses, yet the rationale for this cutoff\nremains insufficiently justified. While the authors briefly explored sensitivity to alternative\nthresholds (8% and 12%), these appear arbitrarily selected and do not constitute a\nsystematic evaluation. Moreover, Figure 2C shows that the freshwater content distribution\nshifts substantially across years, suggesting that a fixed threshold may capture biologically\ndistinct subsets of individuals at different timepoints. This raises concerns about the\ncomparability and interpretation of “jackpot” status across temporal scales. The authors\nare encouraged to provide a stronger conceptual and statistical justification for the 10%\nthreshold, and to consider whether a relative or quantile-based definition (e.g., top x%)\nmight yield more robust or interpretable results given the temporal dynamics of freshwater\nallele accumulation.\nThe 10% cut off was based on the bimodal distribution of the samples collected in 2014, which\nsuggested to us that two distinct groups were in our sample (Review Figure 6). To make this\nclearer, we have updated the visualization of the distribution of freshwater adaptive alleles in\neach year, which was originally presented as a violin plot in supplementary figure S2 (Review\nFigure 6) with the scatter plot in Review Figure 7, as this shows the distinct groups sampled in\nSC2014 more clearly. We include Review Figure 7 in a multi-panel Extended Figure (see\nbelow).\nHowever, we recognize the reviewer’s strong point that our original approach was a subjective\ngrouping and that may introduce biases that affect our downstream interpretation (although we\nnote that the definition of jackpot versus non-jackpot is simply an interpretive tool and is not\nactually implemented in any downstream analytical frameworks, other than the analysis\ncomparing the AFS of jackpot and non-jackpot individuals for the 2014 time point). Therefore,\nto make this inference more robust we have distinguished these two potential groups using a\nmore explicit model-based approach.\nFirst we performed a dip test to check the modality of the SC2014 sample based on freshwater\ncontent (i.e. number of freshwater alleles per individual) and identified two modes (Review\nFigure 8).\nSecondly, we pooled all samples across timepoints for the same metric of the number of\nfreshwater alleles per individual, and used a Bayesian Gaussian Mixture Modelling (BGMM)\napproach implemented in scikit-learn (version 1.7.1) in python 3.11.13 to estimate the number of\nunderlying distributions that best fit the data. BGMM incorporates a Dirichlet process prior that\nallows the model to infer the effective number of components directly from the data. We set the\nmaximum number of components to four and used the full covariance type, which treats each\ncomponent to have its own general covariance matrix. However, doubling the number of\nmaximum components did not change the result of the analysis. Bayesian information criterion\n(BIC) suggested two components to be the best choice, consistent with the dip test (Review\nFigure 9). The fitted model provided posterior probabilities for each observed individual’s\nmembership in each of the two main components based on their freshwater content.\nWe then categorized all the datapoints at each timepoint into one of these two components\n(Review Figure 10). We obtained a remarkably similar profile of individuals assigned to the two\ncomponents at each time point compared to that estimated using our original 10% cut off\n(Review Figure 11), with a Kolmogorov–Smirnov test suggesting that the two distributions were\nnot statistically different (K-S test: statistic = 0.143, p-value = 1). Thus the two components\nappear to reflect the two categories that we previously called non-jackpot carriers and jackpot\ncarriers using the 10% threshold (Review Figure 11). The model categorized individuals with 5%\nfreshwater adaptive alleles as jackpot carriers with a posterior probability of 0.95. Individuals\nwith 8% freshwater adaptive alleles were classified as jackpot carriers with 100% certainty\n(Review Figure 12).\nWe have now adjusted all our results and analysis using this new BGMM-based definition of\njackpot and non-jackpot individuals, though as it results in assignments very similar to the\noriginal ones, it does not change any of our downstream interpretation. In particular, we have\nupdated Fig. 4 and Fig. S17 to reflect the assignment changes. The major inference from Fig. 4\nthat involved jackpot and non-jackpot classification was the suggestion that jackpot carriers\nmated with non-jackpot individuals early on in the colonization process. After applying the\nBGMM classification, this inference is still valid as we observed relatedness between the two\ngroups. The SFS of jackpot vs non-jackpot carriers in Fig. S17 did not change either. We also\nrepeated our Wright-Fisher simulations using a 5% cut off based on the BGMM model rather\nthan the original 10% cutoff, with our initial result that jackpots are necessary for rapid\nfreshwater adaptation still being robust.\nNote that we have now removed all reference to using 10% cutoff in the manuscript as well as\nour previous testing of alternate thresholds and rely completely on the assignments from the\nBGMM which have a more robust statistical justification, as requested by the reviewer. Overall\nwe thank the reviewer for pushing us to address this issue and believe our inference is on firmer\nground because of it.\nIn addition to these changes, we have included/modified the following paragraphs to the results\nand methods to indicate these additional analyses:\nExtended Figure 1 showing the distribution of freshwater adaptive alleles across the years and\nthe BGMM predictions for each individual.\nResults:\n“Therefore, to aid our understanding of genomic change during adaptation to freshwater, we\ncategorized individuals as jackpot carriers and non-jackpot individuals using Bayesian Gaussian\nMixture Modeling (BGMM) (Methods, Supplementary Section 9, Extended Fig. 1B). The BGMM\ncategorized individuals with 5% freshwater adaptive alleles as jackpot carriers with a posterior\nprobability of 0.95 and those with 8% freshwater adaptive alleles were classified as jackpot\ncarriers with 100% certainty (Fig. S4). There were no jackpot carriers in RS2009 and RS2019\nsamples and one jackpot carrier in SC2013 sample (i.e., 1%). There were 23 out of 48 in SC2014\n(i.e., 48%), 84 out of 96 in SC2015 (i.e., 88%), and all in SC2017 and SC2020 (i.e., 100%, Table\nS1). We observed that 27 out of the 47 individuals in SC2014 had contiguous blocks with\nfreshwater-adaptive alleles. These contiguous blocks were mostly heterozygous (Fig. S6) and\nspanned 0.07Mb [0.044cM] to 22.33Mb [70.48cM], with a mean of 2.40Mb [2.9cM] (Fig. 2A,\nSupplementary table 2).”\nMaterials and methods:\n“Classification of individuals as jackpot carriers: We performed a dip test to check the modality\nof the SC2014 sample based on freshwater content (i.e. number of freshwater alleles per\nindividual) and identified two modes. We pooled all samples across timepoints for the same\nmetric of the number of freshwater alleles per individual, and used a Bayesian Gaussian Mixture\nModelling (BGMM) approach implemented in scikit-learn (version 1.7.1) in python 3.11.13 to\nestimate the number of underlying distributions that best fit the data. BGMM incorporates a\nDirichlet process prior that allows the model to infer the effective number of components directly\nfrom the data. We set the maximum number of components to four and used the full covariance\ntype, which treats each component to have its own general covariance matrix. However,\ndoubling the number of maximum components did not change the result of the analysis. Bayesian\ninformation criterion (BIC) suggested two components to be the best choice, consistent with the\ndip test (Fig. S7). The fitted model provided posterior probabilities for each observed\nindividual’s membership in each of the two main components based on their freshwater content.\nWe then categorized all the datapoints at each timepoint into one of these two components (Fig.\nS7). The model categorized individuals with 5% freshwater adaptive alleles as jackpot carriers\nwith a posterior probability of 0.95. Individuals with 8% freshwater adaptive alleles were\nclassified as jackpot carriers with 100% certainty.”\nReview Figure 6\nReview Figure 7\nReview Figure 8\nReview Figure 9\nReview Figure 10\nReview Figure 11\nReview Figure 12\nExtended Figure 1\nRevised Figure 4\nRevised Fig. S17\nRevised Figure 6\n3. The authors acknowledge that elevated homozygosity resulting from inbreeding may\nhave facilitated the purging of maladaptive alleles, thus potentially accelerating positive\nselection at freshwater-adaptive loci. However, this response does not fully address the\npotential mechanistic implications of increased homozygosity on haploblock structure and\nallele frequency dynamics. For instance, greater homozygosity could reduce effective\nrecombination, leading to the apparent expansion or preservation of extended haplotypes\nover time. Additionally, in small populations, hitchhiking effects may be more pronounced,\nallowing non-adaptive alleles to rise in frequency due to linkage with adaptive loci. These\npossibilities warrant at least brief acknowledgment in the discussion to avoid over-\nsimplification of the observed genomic patterns as purely selection-driven. A more nuanced\nconsideration of demographic and genetic drift effects would strengthen the interpretation\nof haploblock dynamics over time.\nWe believe that we are a little bit clearer about the point the reviewer is making here and\nappreciate the need for increased nuance in our text. However we note that there is an important\ndistinction between the increased rate of inbreeding and the increase in freshwater haploblock\nfrequency. The former is indeed likely resulting in increased homozygosity across the entire\ngenome as measured by increased F values. Recessive deleterious recessive variants may be\npurged as a result of this process genomewide, as has been seen in other species going through\nsuch extreme bottlenecks with inbreeding. However the freshwater-adaptive haploblocks are\nactually increasing such that the adaptive loci are transitioning from being primarily homozygous\nfor the marine haplotype to being heterozygous for both freshwater and marine haplotypes, as\nrepresented by the long stretches of yellow in Figure 1. Only much later do such loci start to\nbecome blue, representing homozygosity of freshwater haploblocks, and even then they are less\ncommon than heterozygous haploblocks. Thus effective recombination rates (i.e. observable\nrecombination events) would not be reduced during this stage, and the rise in freshwater\nhaploblocks must therefore be attributed to differential increases in fitness of freshwater\nhaploblocks (i.e positive selection). It is true however that while recessive deleterious alleles in\nmost of the genome may be purged due to inbreeding, deleterious alleles in and around\nfreshwater haplotypes may hitchhike, counterbalancing to some extent the purging of genetic\nload elsewhere in the genome. We have edited/added the following paragraphs to the discussion\nto reflect these complex dynamics.\n“Previous studies have proposed that population bottlenecks may be important during\nfreshwater adaptation by oceanic Threespine Stickleback1,2. These bottlenecks have been\ninferred from an observed reduction in genetic diversity in established freshwater populations1–3.\nHere, we present the first direct observation of the temporal dynamics of stickleback adaptation\nto freshwater during the first nine years of colonization. We observed a bottleneck that occurred\nduring the third and fourth year after founding, likely as a result of reduced fitness of non-\njackpot individuals that made up the majority of the founding population. The resulting reduction\nin population size coincided with a major shift in the genetic composition of the individuals that\nremained in the lake, indicating the population decline was non-random and driven by strong\npositive selection for a small number of jackpot carriers.\nHowever, population growth from a few related individuals will decrease genetic variation and\npotentially increase the population’s genetic load. Genetic load can either be expressed\nphenotypically or masked4 in diploid populations, while bottlenecks can increase phenotypic\nexpression of a population’s genetic load4–6 by exposing rare, recessive deleterious alleles\nnormally masked in heterozygotes. In addition, inbreeding after the bottleneck could increase the\nfrequency of rare deleterious alleles and their probability of homozygosity4,5,7. Such a scenario\ncan lead to one of two outcomes: a) the population could be at risk of extinction due to the\nincrease in the fixation probability of recessive deleterious mutations8–10 through drift, or b)\nnegative selection could be more effective at eliminating slightly deleterious mutations, as\nselection is more potent against recessive or partially recessive mutants in the homozygous\nstate11.\nIn Scout Lake, the latter outcome could, therefore, increase the population’s overall adaptability\nby removing deleterious variants that circulated in the anadromous ancestor alongside the\npositive selection of freshwater-adaptive alleles, at least during the early stages of the adaptive\nprocess, before new deleterious mutations emerge. Such a process has been observed in many\nother species undergoing extreme bottlenecks12–14. Although inbreeding may have led to a\ngeneral purging of deleterious mutations genome-wide, it is interesting to contrast its effects at\nadaptive loci with regards to the dynamics of freshwater haploblocks. While much of the genome\nis becoming increasingly homozygous during these earliest stages of freshwater adaptation,\nadaptive loci are becoming more genetically diverse as homozygous oceanic haploblocks are\ntransitioning into heterozygous haploblocks as the freshwater alleles are increasingly positively\nselected for. This may allow deleterious alleles found within or closely linked to freshwater\nhaplotypes to hitchhike to higher frequencies, thus counteracting the decrease in load elsewhere\nin the genome. The low recombination rates keeping such haploblocks intact over long periods\nmay act somewhat like inversions, which indeed have been shown to accumulate genetic load in\nother species such as Heliconus15. This may even slow the speed of fixation of freshwater loci\nlater in the adaptive process as these deleterious alleles become homozygous, again analogous\nto frequency-dependent selection observed in Heliconus, though unlike inversions, the presence\nof recombination hotspots between adaptive loci may allow some rescue from high genetic load.\nIndeed, even in established populations we rarely see individuals who are fully homozygous for\nall freshwater alleles and there are now emerging empirical studies that explore the dynamics of\ndeleterious alleles linked to adaptive loci in low recombination regions (such as in threespine\nsticklebacks)16,17.”\n4. While the authors reasonably argue that increasing homozygosity could obscure\nrecombination signals, the potential erosion of large freshwater haploblocks remains an\nimportant consideration. Even in the absence of parent-offspring trios, population-level\npatterns, such as reduced haploblock lengths, internal heterogeneity, or linkage\ndisequilibrium decay, could provide indirect but informative evidence of recombination.\nCurrently, the manuscript presents no empirical analysis to support the conclusion that\nrecombination played a limited role. Given that recombination is a central process\ninfluencing haplotype integrity and adaptive allele assembly, a lack of data-driven support\nweakens the strength of this claim. A brief analysis or at least an expanded discussion of\npotential haploblock erosion would strengthen the interpretation that haploblock\npersistence reflects biological processes (e.g., selection, demographic dynamics) rather than\nmethodological limitations or overlooked recombination events.\nWe have now performed statistical tests to determine the differences in overall haploblock\nlengths from one timepoint to another. We first used Welch’s t-test to determine whether the\nmean haploblock lengths are significantly different from one timepoint to another. There was no\nsignificant difference in the mean haploblock lengths between the two Rabbit Slough samples\n(RS2009 vs RS2019; t = 1.2021, p-value = 2.3111e-01), as well as between the RS2019 and\nSC2013(t= 0.4335, p-value = 6.6493e-01). However, there was a statistically significant\ndifference between the mean haploblock lengths of SC2013 and SC2014 (t= -9.7188, p-value =\n3.4958e-21), SC2014 and SC2015 samples (t= -2.2363, p-value = 2.5478e-02) and SC2015 and\nSC2017 (t = -2.7869, p-value = 5.3432e-03). There was no statistically significant difference\nbetween the mean haploblock lengths of SC2017 and SC2020 (t = -1.4509, p-value = 1.4724e-\n01). We further explored the differences in haploblock lengths using Cohen’s D, which expresses\nthe mean difference relative to the pooled standard deviation. We estimated Cohen’s D and\nbootstrapped confidence interval for each pair of consecutive timepoints. There was a small\ndifference between RS2009 and RS2019 (D=-0.1585[-0.352 – 0.087] ) and RS2019 and SC2013\n(D= 0.0044[-0.1341 –0.26] ) compared to SC2013 and SC2014 (D=0.4025[0.367 – 0.4510]).\nBetween SC2014 and SC2015, there was a small difference in the mean haploblock lengths\n(D=0.093[0.012 – 0.172]) compared to SC2013 and SC2014, suggesting that despite the\nstatistically significant t-test, the difference was not substantial at this time period. We observed\nsimilar Cohen’s D between haploblock means in subsequent years: SC2015 and SC2017\n(D=0.083[0.024 – 0.138]); SC2017 and SC2020 (D=0.074[-0.025 – 0.174]). These results\nsuggest that there was limited increase in the haploblock lengths across the time periods as\nmediated by recombination, especially after the bottleneck in 2014.\nAs suggested by the reviewer we have also estimated linkage disequilibrium (LD) decay in and\nacross the freshwater adaptive loci. As previously described for other analyses, we focus on the\n280 loci without any missing data. For each time point, we estimated pairwise r2 among all\ncombinations of the 280 loci using hard-called locus-based genotypes, coded as 0 for\nhomozygous freshwater, 1 for heterozygous, and 2 for homozygous marine, corresponding to the\nnumber of marine alleles carried by each individual at each locus. To quantify the uncertainty in\nLD estimates, we performed bootstrapping using genotype matrices for each of the timepoints.\nFor each timepoint, the genotype matrix was resampled with replacement 1,000 times to generate\na distribution of r2 values in order to estimate 95% confidence intervals (Review Figure 13,\n14\nand 15). We also estimated the r2 between individual SNPs (rather than using the genotype call\nfrom the whole adaptive region) found in the 344 adaptive regions based on the imputed\ngenotypes, and plotted the mean and standard error (Review Figure 16, 17 and 18). We then\nbinned pairs of loci/SNPs by distance to examine differences in the rate of LD decay.\nAs RS2009 and SC2020 are much higher coverage than the other population samples (while also\ncontaining less samples), they will have more accurate genotype calls compared to the imputed\nlow coverage samples, particularly at true heterozygote sites, leading to systematic difference in\nr2 that reflect genotype call accuracy. This effect is visible when comparing the LD decay curves\nbased on individual SNPs for RS2009 vs RS2019, which shows a slight elevation for the former\n(Review Figure 18). Therefore for the individual SNP-based LD analysis RS2009 and SC2020\nwere analyzed separately. However, the loci-based estimates of LD will be more robust to such\neffects as they utilize our approximate likelihood calling approach that leverages information\nfrom multiple SNPs at each loci, and indeed RS2009 and RS2019 largely overlap, and thus both\nlow and high coverage genomes are analyzed together in this case.\nWe have also attempted to test for evidence of partial block loss (haploblock erosion) in our data\nusing an approach that tried to leverage the expected number of meiosis separating biologically\nrelated individuals in our dataset. Unfortunately, due to the low coverage nature of our data, it\nwas not possible to distinguish true haplotype disintegration from individual genotype\nuncertainty.\nWe have updated the text with the following sentences to reflect these observations:\nWe also include Review Figure 15 as an extended figure.\nResults:\n“There was no significant difference in the mean haploblock lengths between the two Rabbit\nSlough samples (RS2009 vs RS2019; t = 1.2021, p-value = 2.3111e-01), as well as between the\nRS2019 and SC2013 (t= 0.4335, p-value = 6.6493e-01). There was a statistically significant\ndifference between the mean haploblock lengths of SC2013 and SC2014 (t= -9.7188, p-value =\n3.4958e-21), SC2014 and SC2015 samples (t= -2.2363, p-value = 2.5478e-02) and SC2015 and\nSC2017(t = -2.7869, p-value = 5.3432e-03). There was no statistically significant difference\nbetween the mean haploblock lengths of SC2017 and SC2020 (t = -1.4509, p-value = 1.4724e-\n01). We further explored the differences in haploblock lengths using Cohen’s D, which expresses\nthe mean difference relative to the pooled standard deviation. We estimated Cohen’s D and\nbootstrapped confidence interval for each pair of consecutive timepoints. There was a small\ndifference between RS2009 and RS2019 (D=-0.1585[-0.352 – 0.087] ) and RS2019 and SC2013\n(D= 0.0044[-0.1341 –0.26] ) compared to SC2013 and SC2014 (D=0.4025[0.367 – 0.4510] ).\nBetween SC2014 and SC2015, there was a small difference in the mean haploblock lengths\n(D=0.093[0.012 – 0.172]), suggesting that despite the statistically significant t-test, the\ndifference was not substantial at this time period. We observed similar Cohen’s D between mean\nhaploblock means between SC2015 and SC2017 (D=0.083[0.024 – 0.138]) as well as SC2017\nand SC2020 (D=0.074[-0.025 – 0.174]). These results suggest that there was limited increase in\nthe haploblock lengths across the time periods, especially after the bottleneck in 2014, which\nwould be expected if the haploblocks were being reassembled by recombination.\nVenu et al. 18 indicated that as many as half of all chromosomes in male, as well as a third of\nchromosomes in females, could be inherited without any recombination event per generation.\nIn addition, recombination is suppressed in individuals that are heterozygous for marine and\nfreshwater haplotypes at adaptive loci. Our sampling period covered nine years, spanning 5-6\ngenerations, and the majority of the individuals sampled after the bottleneck were predominantly\nheterozygous at freshwater adaptive loci. We estimated linkage disequilibrium (LD) across\nfreshwater adaptive loci (Methods) to assess whether recombination had any substantial effect\non the integrity of haploblocks. We observed an elevated levels of LD from SC2014 decay\ncompared to SC2013 and samples from the ancestral population, consistent with a strong\ndirectional positive selection increasing the frequency of adaptive haploblocks or potentially the\neffect of demographic bottleneck that occurred in 2014 (Extended Fig. 2). In subsequent years\nafter the bottleneck, we found no observable patterns of reduction in the rate of LD decay,\nsuggesting that there was not enough time for recombination to have substantial effect on the\nintegrity of haploblocks.”\nMethods\n“Linkage disequilibrium (LD) decay: To estimate the LD decay across the sampling period, we\nfocused on the 280 adaptive loci without any missing data. For each time point, we estimated\npairwise r2 among all combinations of the 280 loci using hard-called locus-based genotypes,\ncoded as 0 for homozygous freshwater, 1 for heterozygous, and 2 for homozygous marine,\ncorresponding to the number of marine alleles carried by each individual at each locus. To\nquantify the uncertainty in LD estimates, we performed bootstrapping using genotype matrices\nfor each of the timepoints. For each timepoint, the genotype matrix was resampled with\nreplacement 1,000 times to generate a distribution of r2 values in order to estimate 95%\nconfidence intervals. We then binned pairs of loci/SNPs by distance to examine differences in the\nrate of LD decay.”\nReview Figure 13\nReview Figure 14\nReview Figure 15\nReview Figure 16\nReview Figure 17\nReview Figure 18\nReferences\n1. Bell, M. Evolution of phenotypic diversity in Gasterosteus aculeatus superspecies on the\nPacific coast of North America. Systematic Biology 25, 211–227 (1976).\n2. Withler, R. E. & Mcphail, J. D. Genetic variability in freshwater and anaromous\nsticklebacks (Gasterosteus aculeatus) of southern British Columbia. Can J. Zool 63, 528–\n533 (1985).\n3. Taylor, E. B. & McPhail, J. D. Historical contingency and ecological determinism interact\nto prime speciation in sticklebacks, Gasterosteus. Proc. Biol. Sci. 267, 2375–2384 (2000).\n4. Bertorelle, G. et al. Genetic load: genomic estimates and applications in non-model animals.\nNat. Rev. Genet. 23, 492–503 (2022).\n5. Crow, J. F. Genetic Loads and the Cost of Natural Selection. in Mathematical Topics in\nPopulation Genetics (ed. Kojima, K.-I.) 128–177 (Springer Berlin Heidelberg, Berlin,\nHeidelberg, 1970).\n6. Wang, J., Hill, W. G., Charlesworth, D. & Charlesworth, B. Dynamics of inbreeding\ndepression due to deleterious mutations in small populations: mutation parameters and\ninbreeding rate. Genet. Res. 74, 165–178 (1999).\n7. Nietlisbach, P., Muff, S., Reid, J. M., Whitlock, M. C. & Keller, L. F. Nonequivalent lethal\nequivalents: Models and inbreeding metrics for unbiased estimation of inbreeding load.\nEvol. Appl. 12, 266–279 (2019).\n8. Hedrick, P. W. Purging inbreeding depression and the probability of extinction: full-sib\nmating. Heredity 73 ( Pt 4), 363–372 (1994).\n9. Lynch, M., Conery, J. & Burger, R. Mutation Accumulation and the Extinction of Small\nPopulations. Am. Nat. 146, 489–518 (1995).\n10. Palkopoulou, E. et al. Complete genomes reveal signatures of demographic and genetic\ndeclines in the woolly mammoth. Curr. Biol. 25, 1395–1400 (2015).\n11. Barrett, S. C. H. & Charlesworth, D. Effects of a change in the level of inbreeding on the\ngenetic load. Nature 352, 522–524 (1991).\n12. Grossen, C., Guillaume, F., Keller, L. F. & Croll, D. Purging of highly deleterious\nmutations through severe bottlenecks in Alpine ibex. Nat. Commun. 11, 1001 (2020).\n13. Croft, L., Matheson, P., Butterworth, N. J. & McGaughran, A. Fitness consequences of\npopulation bottlenecks in an invasive blowfly. Mol. Ecol. e17492 (2024).\n14. Dussex, N., Morales, H. E., Grossen, C., Dalén, L. & van Oosterhout, C. Purging and\naccumulation of genetic load in conservation. Trends Ecol. Evol. 38, 961–969 (2023).\n15. Jay, P. et al. Publisher Correction: Mutation load at a mimicry supergene sheds new light on\nthe evolution of inversion polymorphisms. Nat. Genet. 53, 763 (2021).\n16. Harringmeyer, O. S. & Hoekstra, H. E. Chromosomal inversion polymorphisms shape the\ngenomic landscape of deer mice. Nat. Ecol. Evol. 6, 1965–1979 (2022).\n17. Nickel, J., Laine, J. & Foote, A. D. Complex patterns of hitchhiking mutation load among\nstickleback populations. bioRxiv (2025) doi:10.1101/2025.07.24.666323.\n18. Venu, V. et al. Fine-scale contemporary recombination variation and its fitness\nconsequences in adaptively diverging stickleback fish. Nat. Ecol. Evol. 8, 1337–1352\n(2024)."}, {"round": 3, "reviewer_comments": {"reviewer_1": "The authors have fully addressed my concerns by adding multiple robust new analyses and substantially strengthening the methodological and mechanistic interpretation. The revisions are thorough and well executed, and I have no further comments."}, "author_response": ""}]

Decision letter