The genomics of postglacial vicariance and freshwater adaptations in European subarctic threespine sticklebacks

Marine threespine sticklebacks (Gasterosteus aculeatus) have colonized numerous freshwater lakes since the last ice age. The loss of body armor, such as reduced pelvic spines and lateral plate numbers, is a recurrent feature upon freshwater colonization and is attributed to parallel evolution. This study examines genome-wide genetic diversity and differentiation among sticklebacks from a marine site and two freshwater lakes of the same watercourse in subarctic Europe. The upper lake is unique in that it harbors a population of polymorphic sticklebacks, some with fully developed and others with reduced pelvic structures. Our results based on deep sequencing of pooled population samples showed common signatures of selection for freshwater sticklebacks in certain parts of the genome, such as the Eda containing region, but also evidence of differential selection, and the presence of large chromosomal inversions that seem to play an essential role in stickleback evolution. Pelvic reduction in sticklebacks has previously been linked to deletions in the enhancers of the pituitary homeobox transcription factor gene (Pitx1). While the genetics of Pitx1 seem unable to fully explain pelvic spine polymorphism in this population, we found differentiation between spined and spineless sticklebacks in several genomic regions, which harbor genes that might be involved in pelvic development. Most significantly, genetic differentiation between spined and spineless sticklebacks was noted in a region of chromosome 9 where the gene Hand2, previously implicated in limb development, is located. Our findings suggest that pelvic reduction in these sticklebacks involves multiple genetic factors, indicating parallel evolution through polygenic influences.

1 Introduction

The threespine stickleback (Gasterosteus aculeatus) is a small teleost fish species (typically 30-60 mm in length) that inhabits coastal and inland waters in the Northern Hemisphere. Marine sticklebacks are anadromous; they live in marine habitats, but migrate to breed in brackish or fresh waters. Being tolerant to changes in salinity, marine threespine sticklebacks commonly colonize freshwater bodies and establish novel populations, some of which, upon restricted gene flow, form freshwater ecotypes that complete their entire life cycle in freshwater. The evolutionary trajectory from marine to freshwater ecotypes in threespine sticklebacks involves a range of changes in life history, physiology, behavior, and morphology, among which the reduction of external body armor is a common and recurrent feature. Marine sticklebacks invariably possess strong armor, including lateral bony plates, three dorsal spines, and two pelvic spines, which serve as protection against gape-limited predators such as birds and fishes (Wootton, 1976). In contrast, freshwater sticklebacks typically exhibit fewer bony plates and, less commonly, reduced pelvic structures (Figure 1). These reductions in body armor are thought to be due to changes in selective regime that sticklebacks experience upon transitioning to freshwater habitats, which often feature lower calcium availability and relaxed predation pressure compared to marine environments (Bell et al., 1993; Hagen and Gilbertson, 1973; Klepaker, 1995; Moodie and Reimchen, 1976; Reimchen, 1983).

Figure 1

Figure 1. Variation in morphological features in threespine sticklebacks. (A) A photographic image of a marine threespine stickleback showing dorsal spines and pelvic spines. Schematic diagrams of (B) marine stickleback with well-developed lateral plates, dorsal spines, and pelvic spines; (C) freshwater stickleback with reduced lateral plates and well developed pelvic spines; (D) freshwater stickleback with a complete reduction of pelvic spines and reduced lateral plates.

In principle, the evolution of adaptive traits in organisms colonizing new habitats could arise either from selection on standing genetic variation (SGV) in the ancestral population or from de novo mutations in the founder population (reviewed by Schluter and Conte, 2009). SGV refers to the pre-existing genetic variation, or the presence of multiple alleles at particular loci, within a population. Several authors have argued that alleles from SGV are more commonly involved than de novo mutations for several reasons: (i) they are usually present in higher frequencies, (ii) they are immediately available in the new habitat, and (iii) they have already been tested in past environments (Barrett and Schluter, 2008; Innan and Kim, 2004; Schluter et al., 2004). Supporting this view, there are indications of the reuse of ancestral genetic variants for freshwater adaptation in sticklebacks, distributed across multiple chromosomes and genomic regions (Hohenlohe et al., 2010; Jones et al., 2012; Liu et al., 2018). Some of these genetic variants, including some that reside within chromosomal inversions, appear to be ancient, dating back several million years (Nelson and Cresko, 2018). However, it has also been shown that the evolution of similar phenotypic traits in subspecies or populations of the same species could arise either from de novo mutations in genes known to be involved, or from the combined effects of several loci, as the traits in question are often quantitative and polygenic (reviewed by Arendt and Reznick, 2008; Hoekstra et al., 2006).

The number of ancestral alleles within the SGV of sticklebacks that are adapted to freshwater environments could be reduced due to genetic drift and selection against these alleles during range expansions through the marine environment. Thus, we expect the amount of shared genetic variants for freshwater adaptation to be associated with colonization history and geographic distance among a given collection of populations. Fang et al. (2018) reconstructed the worldwide phylogenetic relationships and colonization history of threespine sticklebacks, inferring that current populations originated in the Pacific Ocean during the late Pleistocene. Sticklebacks then colonized the Atlantic Ocean through the Arctic approximately 40-50 thousand years ago, initially forming a southern European clade and later a derived trans-Atlantic clade that includes stickleback populations in eastern North America and northern Europe. In compliance with this, shared ancestral polymorphisms appear to be more common among Eastern Pacific freshwater locations than on a global scale (DeFaveri et al., 2011; Fang et al., 2020; Jones et al., 2012). Also, while the majority of earlier studies were based on limited geographic sampling and focused on the Eastern Pacific region, recent studies have suggested the presence of previously unexplored large-scale geographic heterogeneity in the genomic basis of parallel evolution among sticklebacks (Fang et al., 2020; Terekhanova et al., 2019).

The loss of lateral bony plates, a feature most typical of marine-freshwater transitions in sticklebacks, can evolve within a few decades (Bell, 1994; Klepaker, 1993; Roberts Kingman et al., 2021). The ectodysplasin gene (Eda) on chromosome 4, which is regulated by a cis-regulatory element, appears to be responsible for the majority of variation in bony plates (Colosimo et al., 2004; Cresko et al., 2004; O’Brown et al., 2015). However, several loci of minor importance, mapping to other linkage groups, have also been implicated (Colosimo et al., 2004; Cresko et al., 2004; Peichel et al., 2001). The reduction of lateral bony plates in freshwater populations is commonly attributed to selection on SGV (Colosimo et al., 2005; Schluter and Conte, 2009).

While the loss of lateral bony plates is common in freshwater populations, the partial or complete loss of pelvic spines is less frequent, except in certain freshwater populations such as those in Cook Inlet, Alaska (Bell et al., 1993; Bell and Ortí, 1994; reviewed by Klepaker et al., 2013). A few other populations with reduced pelvic structures are found in lakes across British Columbia, Iceland, Scotland, and Norway (Chan et al., 2010; Coyle et al., 2007; McPhail, 1992; Peichel et al., 2001; Shapiro et al., 2004; reviewed by Klepaker et al., 2013). In Norway, pelvic reduction is reported in only eight out of more than 200 examined populations (Klepaker and Østbye, 2008; Klepaker et al., 2013). A major determinant of pelvic development in threespine sticklebacks is the pituitary homeobox transcription factor gene Pitx1, located on chromosome 7 (Chan et al., 2010; Coyle et al., 2007; Cresko et al., 2004; Shapiro et al., 2004). The upstream enhancer PelA of the Pitx1 is crucial for regulating its expression, and the presence of multiple TG-repeats in PelA makes it prone to mutation, providing a molecular mechanism for recurrent pelvic reductions (Chan et al., 2010). Another enhancer of Pitx1, the PelB, and additional loci on chromosomes 2, 4, and 8 are also suggested to be involved in the development of pelvic spines (Peichel et al., 2001; Shapiro et al., 2004; Thompson et al., 2018).

In this study, we examined the genomics of polymorphic sticklebacks from two lakes within the same watercourse, along with specimens from a nearby marine location in subarctic Norway, Northern Europe. These sticklebacks are presumably among the most distantly related to ancestral Pacific populations (Fang et al., 2018). Sticklebacks from both lakes exhibit the typical freshwater phenotypic feature of having fewer lateral plates than marine conspecifics. Additionally, sticklebacks from the upper lake of these two lakes exhibit polymorphism in their pelvic structures, ranging from fully spined to asymmetrically spined and spineless. Previously, we showed that the number of TG-repeats in PelA differed between the two lakes with PelA being consistently shorter in individuals from the upper lake. However, no clear association was observed between PelA enhancer variants and pelvic status among individuals from the upper lake. These findings suggest that additional loci may be involved in the variable manifestation of pelvic spines in this species (Adhikari et al., 2023).

We sequenced pooled DNA samples of sticklebacks from each of the three study locations to provide a first assessment of genome-wide diversity and differentiation, including chromosomal inversions, among marine and vicariant freshwater sticklebacks in the European subarctic. We sought to identify genomic regions displaying signatures of selection potentially involved in the adaptation of sticklebacks to freshwater environments. Further, we conducted a bulk segregant analysis of spined and spineless sticklebacks from the same lake, specifically aiming to identify differences in the Pitx1 locus and genomic signatures of differential selection between the two groups, and to indicate candidate genes that could be involved in the evolution of pelvic reduction. Two individual sticklebacks, one spined and one spineless specimen, were further subjected to long-read sequencing, aiming to resolve their Pitx1 genotypes.

2 Materials and methods

2.1 Sampling and data collection

A total of 426 threespine sticklebacks were collected from two freshwater lakes - Lake Storvatnet (68°46′49″N, 15°9′36″E; 80 m altitude), and Lake Gjerdhaugvatnet (68° 46’17’’N, 15° 9’2’’E; 20 m altitude), which are part of the same watercourse (Figure 2) in June 2017, 2019 and 2020. Additional samples were taken from a nearby marine location at the tidal mouth of a small river at Sandstrand (68°44′45″N, 15°20′42″E), Langøya Island in Northern Norway in June 2020 (Figure 2). Several waterfalls between the two lakes and between the lower Lake Gjerdhaugvatnet and the sea, prevent gene flow between the three stickleback populations. The sticklebacks from the marine site were collected at the outlet of a small river, where they can move freely between water of varying salinity, which is dynamic due to the tidal cycle. The total length of each specimen was measured, and only specimens ≥ 3 cm long were included. The specimens were euthanized with tricaine methanesulfonate (MS222), and care was taken to preserve their pelvic structures. The caudal fin was cut off and discarded, and 5 mm posterior fin-muscle samples were collected, homogenized by bead beating using a Dremel 8220 rotary tool (MP Biomedicals) and 0.5 ml DNA/RNA Shield solution (Zymo Research), and then kept at low temperatures for subsequent analysis. Morphological traits, such as lateral plate number and pelvic score (PS), were assessed using a 20x magnifying stereomicroscope. Each side of the fish (left and right) was assigned a PS between 0 to 4 (as defined by Bell’s five-graded scale; Bell, 1987). The combined pelvic scores (CPS), ranging from 0 to 8, were then computed (Adhikari et al., 2023) (see Supplementary S1). Based on lateral plate morphology, specimens were categorized as “completely plated”, “partially plated”, and “low plated” (Klepaker et al., 2012; O’Brown et al., 2015; Schröder et al., 2023; Wootton, 1976) (see Supplementary S2).

Figure 2

Figure 2. Sampling sites in Norway. Map showing the sampling sites at Langøya island: the upper Lake Storvatnet (68° 46’49’’N, 15° 9’36’’E), the lower Lake Gjerdhaugvatnet (68° 46’17’’N, 15° 9’2’’E) in the watercourse, and the marine sampling site (68° 44’45’’N, 15° 20’42’’E) located about 8 km from the two lakes.

2.2 Preparation of samples for pooled DNA sequencing

DNA was extracted from 20 µl of stored homogenized samples with a Monarch genomic DNA purification kit (New England Biolabs), and quality checked with Qubit 4 fluorometer (Invitrogen, Thermo Fisher Scientific). DNA samples with concentrations ≥20 ng/µl were included.

Four DNA pools were prepared from three locations, with each pool containing 40 individuals (Karlsen et al., 2024). The four DNA pools were: i) “Spined-Storvatnet”, featuring specimens from Lake Storvatnet with fully developed symmetric pelvic spines on both sides (CPS = 4 + 4), ii) “Spineless-Storvatnet”, also from Lake Storvatnet, including specimens with reduced pelvic spines or pelvic apparatus on both sides (spineless – defined by CPS= 2-4; see Table 1), iii) “Gjerdhaugvatnet”, and iv) “Marine”. Thus, although symmetric-spined, spineless, and asymmetric-spined sticklebacks were sampled in Lake Storvatnet, only DNA pools consisting of spined and spineless sticklebacks were prepared for genomic comparison of the two morphotypes through bulk segregant analysis. No sticklebacks with reduced pelvic spines were found in Lake Gjerdhaugvatnet or the marine site; all specimens from these locations were fully spined (CPS = 4 + 4) (Adhikari et al., 2023).

Table 1

Table 1. Morphological measurements of threespine sticklebacks from the three examined populations.

First, DNA samples with concentrations above 20 ng/µl were diluted to 20 ng/µl. Then 5 µl (100 ng) aliquots from each of 40 DNA samples were mixed to prepare a 200 µl (4000 ng) DNA pool in a 1.5 ml Eppendorf tube (Sham et al., 2002), using the elution buffer provided with the kit as a diluent. The concentration of each DNA pool was verified using both Nanodrop and a Qubit 4 fluorometer, and sample quality was assessed by agarose gel electrophoresis and visual inspection.

2.3 Sequencing

2.3.1 Library preparation and sequencing

The four DNA pools were sent to the Norwegian Sequencing Center (NSC) at the University of Oslo (UiO) for whole genome sequencing (WGS). From each of the four DNA pools, 250 ng genomic DNA was used as input to the Illumina DNA Prep protocol. This protocol applies the tagmentation procedure to fragment DNA and add adapters (formerly known as Nextera Flex; Illumina, San Diego, CA, USA). Libraries were prepared according to the manufacturer’s procedures, using six cycles of PCR amplification and incorporating unique-dual-indexes. Sequencing was performed on a NovaSeq instrument (Illumina) using a ¼ S4 flow-cell, with the XP workflow, employing 150 bp paired-end reads following the manufacturer’s recommendations and running RTA v3.4.4. Demultiplexing was carried out using bcl2fastq v2.20.0.422, and the data were received as raw reads from NSC.

2.3.2 Bioinformatic analysis

Before further processing, the raw data for the pooled population samples with higher read numbers were downsized to approximately the same number of reads as the “Spined-Storvatnet” pool using the reformat tool of BBMap (version 39.01) (https://sourceforge.net/projects/bbmap/). Raw reads were trimmed from adapters and low-quality bases using the Fastp program (version 0.20.0) and data quality was checked using the FastQC program (version 0.11.9). Sequence alignment was performed with the BWA-mem alignment algorithm (version 0.7.17-r1188), producing SAM mapping files. During alignment, the data were mapped to the recent version of the reference genome of threespine stickleback, GAculeatus_UGA_version5 (GeneBank assembly accession: GCA_016920845.1; RefSeq assembly accession: GCF_01692084). The SAM files were converted into sorted BAM files using Samtools version 1.13. PCR duplicates from reads were removed with Picard.jar (version 2.27.2). The statistics, coverage, and depth of the BAM files were checked with Samtools. Subsequently, Mpileup and Synchronization (sync) files were generated from the BAM files, and these files were used with the PoPoolation tool (version 1.2.2 and version 2_1201) (Kofler et al., 2011a, b) to analyze population genetics parameters for the stickleback population samples (see Supplementary S3).

2.4 Estimation of population genetic parameters among freshwater and marine sticklebacks

We used the PoPoolation toolbox (version 1.2.2) (Kofler et al., 2011a) to estimate (i) nucleotide diversity (π), which represents the average number of nucleotide differences per site between DNA sequences in a population sample (Nei and Li, 1979), and (ii) Tajima’s D (T_D), a test statistic for assessing neutrality by comparing the mean pairwise difference (π) between sequences in a population sample to the number of polymorphic sites (s) (Tajima, 1989) across all chromosomes for the three population samples. It should be noted that the “Spined-Storvatnet” pool was used to represent stickleback specimens from Lake Storvatnet, as the other two populations contained only spined sticklebacks. We employed a sliding window approach with a window size of 4000 bp and a step size of 2000 bp (see Supplementary S3; Kofler et al., 2011a).

Genetic differentiation, quantified using F_ST, measures the genetic variance between two or more populations (Gregorius, 1987) based on allele frequencies. In this study, pairwise F_ST values were calculated between the three population samples using the PoPoolation2 tool (version 2_1201) (Kofler et al., 2011b). F_ST values were computed across all chromosomes using a sliding window approach with a window size of 4000 bp and a step size of 2000 bp, consistent with the parameters used for analyzing nucleotide diversity and T_D (see Supplementary S3).

The resulting nucleotide diversity, T_D, and F_ST values were then imported into R (version 4.2.0, R studio v 1.4.1717) for downstream analysis. Our aim was to identify genomic regions displaying structural variation such as inversions, and signatures of selection indicating local adaptation. We also identified loci potentially involved in limb or bone development, residing within genomic regions that were highly differentiated between sites. We specifically focused on examining key loci, including the Eda locus and the Pitx1 locus, to understand their role in the variation of bony armor structures observed across the studied stickleback populations.

2.5 Bulk segregant analysis of spined and spineless sticklebacks from Lake Storvatnet

2.5.1 Whole genome analysis

The genetic basis for pelvic spine reduction in sticklebacks from Lake Storvatnet is not well understood (Adhikari et al., 2023). We calculated nucleotide diversity and T_D along the genome for spined and spineless groups using “Spined-Storvatnet” and “Spineless-Storvatnet” DNA pools, respectively, as previously described. To identify genomic loci associated with pelvic spine reduction, we examined genome-wide differentiation (F_ST) between the spined and spineless groups, applying the same parameters as outlined above. Fisher’s Exact test (with the same parameters as described above) was then employed to estimate statistically significant differences in SNP frequencies between the two groups, with a significance level set to -log₁₀(p) < 5x10^-8 (Kofler et al., 2011b) (see Supplementary S3).

2.5.2 Pitx1 (BAC clone) analysis

Unfortunately, the latest reference genome assembly against which we mapped the PoolSeq reads does not include the Pitx1 locus. Previous reference genome assemblies also lack this locus. The Pitx1 gene is located within a repetitive region at the sub-telomeric end of chromosome 7, which poses challenges for assembling it as a linked contig. Much of the work on the complete sequence of this region is based on Sanger-sequenced BAC libraries (GenBank: GU130435.1) (see Chan et al., 2010). Therefore, to examine the Pitx1 locus, we mapped the PoolSeq data from spined and spineless specimens from Lake Storvatnet against the BAC sequences. As controls, we used PoolSeq data from downstream Lake Gjerdhaugvatnet and the marine site. Coverage depths at the Pitx1 locus were calculated from the mapped BAM files for both spined and spineless sticklebacks from Lake Storvatnet, as well as the samples from Lake Gjerdhaugvatnet and the marine site. Following this, we estimated nucleotide diversity, T_D, and genetic differentiation using the previously described parameters.

2.5.3 Nanopore long-read sequence analysis of the Pitx1 locus

We performed Nanopore long-read sequencing of one spined (CPS=4 + 4 = 8) and one spineless stickleback (CPS=1 + 1 = 2) from Lake Storvatnet to further inspect the Pitx1 locus. The Oxford Nanopore Technologies (ONT) library was prepared using the Native Barcoding Kit 24 V14 (SQK-NBD114.24), following the manufacturer’s standard protocol. Total DNA for ONT sequencing was extracted using the Monarch HMW DNA Extraction Kit (New England Biolabs). Duplex sequencing was performed on the PromethION P2 Solo platform, following the recommended procedures for priming and loading the flow cell.

The raw sequencing data were processed using the Dorado basecaller (v. 0.5.0) to convert raw reads into unaligned BAM format. These unaligned BAM files were then converted into FASTQ format using Samtools (version 1.13). The Porechop (version 0.2.4) was used to trim adapters and chimeric reads present in the FASTQ sequences, and NanoFilt (version 2.6.0) was used to remove low quality reads. The trimmed FASTQ files were mapped with the BAC libraries (GenBank GU130435.1) to inspect changes within Pitx1 locus with the help of the minimap2 (version 2.24-r1122) program to create SAM files. SAM files were then converted to sorted BAM files with the help of Samtools. Further, coverage depth of the BAM files was compared between the two individual (spined S01 and spineless S27) sticklebacks to identify differences along the Pitx1 locus.

2.6 Animal welfare

The study was carried out according to ethical guidelines stated by the Norwegian Ministry of Agriculture and Food through the Animal Welfare Act. According to these guidelines, we were not required to, and therefore do not have, a specific approval or approval number.

3 Results

Of the 304 sticklebacks from Lake Storvatnet, 37% were fully spined, 33% had asymmetric spines (mostly left-biased asymmetric), and 30% were spineless (Table 1). In contrast, all specimens from Lake Gjerdhaugvatnet (N=73) and the marine site (N=50) had fully developed right and left pelvic spines, proportionate to their body size (Figure 3). For this study, only fully spined and spineless specimens from Lake Storvatnet were analyzed. All examined sticklebacks in both Lake Storvatnet and Lake Gjerdhaugvatnet were low plated. In contrast, the marine pooled sample consisted of 32 completely plated, 5 partially plated, and 3 low plated specimens (see Supplementary S1, S2).

Figure 3

Figure 3. Scatter plot presenting the distribution of specimens among three locations (Lake Storvatnet, Lake Gjerdhaugvatnet, and the marine site) based on the total length of the body and the length of left and right pelvic spines. In general, specimens from the marine site have longer pelvic spines relative to their body size compared to those from the lakes. Specimens from Lake Storvatnet exhibit the shortest or most reduced spines in proportion to body size.

3.1 Genome coverage and sequencing depth of pooled DNA samples

Overall sequencing depth (raw data) for each DNA pool ranged from 507x to 728x. Approximately 660 million (M) reads per population sample were successfully mapped against the threespine stickleback reference genome. Of these, about 94% were properly paired, yielding an average mapping quality of 35.7 for each population sample (see Supplementary S4). Across the chromosomes, sequencing coverage ranged from 96% to 99%, with a depth of 159x to 225x, except for the Y chromosome (chr). The Y chromosome showed lower coverage depth in all samples compared to other chromosomes and was therefore excluded from the subsequent calculations of population genetic parameters (see Supplementary S5).

3.2 Genome-wide patterns of genetic diversity and differentiation among populations

The assessment of genome-wide diversity and differentiation among sticklebacks from a marine population and two freshwater lakes in Northern Norway revealed demographic patterns and evidence of selection, shedding light on the evolutionary processes shaping these populations. The marine population exhibited an average nucleotide (genetic) diversity (π) of 0.33% (0.0033), while sticklebacks from Lake Gjerdhaugvatnet and Lake Storvatnet showed lower average π values of 0.21% and 0.23%, respectively (see Supplementary S6, S10). Additionally, the marine population had a more negative average Tajima’s D (T_D) of -1.10, indicating a possible excess of rare alleles (purifying selection), compared to the less negative averages of -0.73 in Lake Gjerdhaugvatnet and -0.78 in Lake Storvatnet (see Supplementary S6, S11). The average genetic differentiation (F_ST) value between sticklebacks from the two freshwater populations was relatively high (0.21), reflecting pronounced genetic separation between these populations. Furthermore, marine sticklebacks exhibited slightly greater genetic differentiation from sticklebacks in Lake Storvatnet (average F_ST = 0.15) compared to those in Lake Gjerdhaugvatnet (average F_ST = 0.14) (Figure 4; see Supplementary S7, S12). Notably, all of the population genetic parameters F_ST, π, and T_D varied widely among chromosomes and chromosomal segments (see Supplementary S6, S7, S10-S12).

Figure 4

Figure 4. Boxplots presenting F_ST values among three populations across 21 chromosomes. F_ST values between Lake Gjerdhaugvatnet and the marine populations, and between Lake Storvatnet and the marine populations, are shown in blue and grey boxes, respectively, while F_ST values between the two freshwater populations (Storvatnet and Gjerdhaugvatnet) are shown in red boxes. Black points within each box indicate the average F_ST values for each chromosome. Notably, the average F_ST values between the two freshwater populations are higher compared to the other two pairs (marine and freshwater specimens) along each chromosome, with the exception of chromosome 21.

3.2.1 Genomic segments showing signatures of selection

Prominent signatures of selection acting on specific parts of several autosomes (chromosomes 1, 3, 4, 5, 6, 7, 8, 9, 11, 12, 14, 17, 20, and 21) were observed based on patterns in π, T_D and F_ST values along the genome (see Supplementary S10-S12). Typically, we observed evidence of directional selection in specific parts of the genome, characterized by highly negative T_D values (~ -2.0 to -4), lower π values, and high genetic differentiation between one or both of the freshwater populations and the marine (exemplified in Figure 5). In specific, on chr 9, regions spanning 12.1-12.3 Mb, and 12.6-12.8 Mb (Figure 5A) show indications of directional selection in Lake Storvatnet. Similarly, in Lake Gjerdhaugvatnet, patterns indicative of directional selection were detected in a region spanning 10.1-10.8 Mb on the same chromosome (Figure 5B). These observations underscore directional selection as an important evolutionary force shaping the genetic landscape of these subarctic freshwater stickleback populations.

Figure 5

Figure 5. Population genetic parameters of selected regions at chromosome 9 indicating signatures of selection. (A) In Lake Storvatnet genomic regions on chr 9 (from 12.1 to12.3 Mb, and 12.6 to 12.8 Mb) characterized by highly negative T_D values (~ -2.0 to -4), lower π values, and higher genetic differentiation compared to other populations. Similarly, (B) in Lake Gjerdhaugvatnet, regions on chr 9 (from 10.1 Mb to 10.8 Mb) are characterized by highly negative T_D values, lower π values, and higher genetic differentiation compared to other populations.

3.2.2 Identification of genes potentially involved in bone and limb development

Based on high genetic differentiation of specific genomic regions between sites, we identified several loci that may be associated with limb or bone development and phenotypic variation. We found highly differentiated regions between the two freshwater populations containing the T-box transcription factor 4 (chr 1; Tbx4; 18.98 to 19.01 Mb), fibroblast growth factor 8a (chr 6; Fgf8a; 27.36 to 27.41 Mb), POU class 1 homeobox 1(chr 16; Pou1f1; 13.13 to 13.14 Mb), and ALX homeobox 1 (chr 19; Alx1; 13.56 to 13.57 Mb). In Lake Storvatnet sticklebacks, regions on chromosome 14, encompassing astrotactin-2-like (Astn2/LOC120831669; 5.32 to 5.48 Mb), bone morphogenetic protein/retinoic acid inducible neural-specific 1 (Brinp1; 5.12 to 5.21 Mb), and LIM homeobox 3 (Lhx3; 4.99 Mb to 5.0 Mb), were found to be highly differentiated compared to the two other sites. Additionally, the 3.0 to 3.4 Mb region on chromosome 5, including the growth/differentiation factor 10-like gene (Gdf10), bone morphogenic protein 2-like gene (LOC120819618/Bmp2), and transcription factor Sox-9-A-like gene (Sox9a), as well as two regions (11.0 to 12.0 Mb and 13.0 to 14.0 Mb) on chromosome 7, containing LOC120822309 (histone H2A.Z: H2az), MLLT3 super elongation complex subunit (Mllt3), LOC120822009 (HMG box transcription factor BBX: Bbx), brain-specific homeobox (Bsx), and LOC120822868 (histone H2A: H2a), exhibited elevated F_ST values between marine and Lake Storvatnet sticklebacks (see Supplementary S14).

3.2.3 Genomic divergence at the Eda locus on chromosome 4

We observed lower genetic diversity (π) across the Eda locus in freshwater populations (Figure 6A). Further, we also observed an elevated F_ST peak at the Eda locus between freshwater and marine specimens, supported by coverage depth assessments of the Eda locus from pooled sequences (Figures 6B, C). Notably, there was a lack of freshwater sequences aligning to specific intronic regions near exon 2, where a possible deletion of ~350 bp in freshwater sticklebacks may contribute to the observed genetic differentiation in this area (Figures 6B, C). Furthermore, we identified nine SNPs within exons 1 and 8 that were fixed for the reference allele in freshwater sticklebacks, while marine sticklebacks were polymorphic, with the reference allele occurring as the minor allele. Similarly, in the intergenic region between Eda and TNFSF13B, three SNPs (Ref/Alt: T/C, G/T, and G/A) were fixed for the reference allele in freshwater sticklebacks, whereas marine sticklebacks were polymorphic (Table 2, Figure 6D).

Figure 6

Figure 6. Illustration of the Eda locus between freshwater and marine sticklebacks. (A) Nucleotide diversity (π) and Tajima’s D plots for the Eda locus. The red line represents the Lake Storvatnet population which shows lower genetic diversity and Tajima’s D values compared to other populations. (B) F_ST plot showing genetic differentiation across the Eda locus on Chromosome 4 among sticklebacks from Lake Gjerdhaugvatnet, Lake Storvatnet, and a marine site. The blue line represents F_ST values for the Gjerdhaugvatnet versus marine comparison, while the red and gray lines represent Storvatnet versus Gjerdhaugvatnet and Storvatnet versus marine, respectively. Higher genetic differentiation was observed between the marine and freshwater sticklebacks at this locus. (C) Coverage depth plot for the same region as shown in A, indicating sequencing depth for each population. The region highlighted in grey in both B and C near to the exon 2 of Eda represents an area with sequence coverage for marine sticklebacks only. (D) Gene structure of the Eda locus, indicating differentiation at a number of nucleotide positions between freshwater and marine sticklebacks. Freshwater sticklebacks were fixed for the reference nucleotide (left) at these positions, while the alternate allele (right) occurred as the major allele in marine sticklebacks (see Table 2).

Table 2

Table 2. Selected regions along the Eda locus (chr 4) showing differentiation at single nucleotide positions between freshwater stickleback populations (Lake Storvatnet and Lake Gjerdhaugvatnet) and marine sticklebacks.

3.2.4 Assessment of genome wide diversity and differentiation suggests chromosomal rearrangements

Our genome-wide analysis revealed potential chromosomal rearrangements, such as inversions, which have been reported in other stickleback populations (Jones et al., 2012). In chr 1, elevated π values between 26.0 to 26.5 Mb were observed in the marine population but were absent in the lake populations (Figure 7A), suggesting a clear differentiation between marine and freshwater sticklebacks. This, combined with high F_ST values between marine and lake populations (see Supplementary S12) and elevated T_D values in the marine sticklebacks (Figure 7B, and see also Supplementary S11), points to either balancing selection or a marine-specific inversion, at least in a subset of marine samples. In chr 11, a region spanning 6.2-6.6 Mb, harboring 23 genes, shows strong peaks in both π and T_D in Lake Storvatnet (Figures 7C, D), indicating the likely presence of both inverted and non-inverted variants, suggestive of balancing selection in this population. Similarly, in chr 21, we identified a potential inversion spanning 9.8-11.7 Mb in both lake populations. The patterns of π (Figure 7E) and T_D (Figure 7F) in this region suggest the presence of both inverted and non-inverted variants, with T_D values in Lake Storvatnet particularly emphasizing the inversion, while F_ST comparisons suggest a more balanced mix of these variants in Lake Storvatnet compared to Lake Gjerdhaugvatnet (Figures 7E–G).

Figure 7

Figure 7. Population genetic parameters of selected regions at chromosomes 1, 11 and 21 indicating putative chromosomal inversions. Chromosome 1 (26 to 26.5 Mb); based on (A) nucleotide diversity pi, and (B) Tajima’s (D) Chromosome 11 (6.2 to 6.6 Mb); based on (C) nucleotide diversity pi, and (D) Tajima’s (D) Chromosome 21 (9.8 to 11.7 Mb); based on (E) nucleotide diversity pi, (F) Tajima’s D, and (G) pairwise F_ST.

3.3 Bulk segregant analysis of spined and spineless sticklebacks from Lake Storvatnet

3.3.1 Whole genome analysis

In our comparative analysis of spined and spineless sticklebacks from Lake Storvatnet, we observed consistent patterns in π and T_D values between the two groups, which is expected given that both groups belong to the same population (see Supplementary S6, S10, S11). The overall genetic differentiation between the groups was low, as reflected by a mean F_ST value of 0.0099, which aligns with the notion of minimal divergence across most of the genome (see Supplementary S7, S13).

We identified 37 SNPs exceeding the genome-wide significance threshold (-log₁₀(p) value 5x10^-8) (Figure 7A). A notable discovery was a highly differentiated region on chromosome 9 (spanning 4.0-4.4 Mb), which we refer to as an “F_ST-island” (Figures 8A, B). This island contains 31 genes (see Supplementary S8), including a gene of particular interest, Hand2, which is known for its role in hind limb development and its potential involvement in pelvic spine formation. Additionally, the spineless group exhibits higher π and elevated T_D values compared to the spined group within the “F_ST-island” (Figure 8C), suggesting that balancing selection may be maintaining heterozygous alleles at higher frequencies in the spineless group.

Figure 8

Figure 8. Manhattan Plots and nucleotide diversity analysis between spined and spineless groups from Lake Storvatnet. (A) Manhattan plot visualizing -log 10 (p) values derived from Fisher’s exact test, identifying significant differences in allele frequencies between the two groups. Thirty-seven SNPs surpass the significance threshold of -log10(p) = 5x10^-8 (blue line), indicating that these loci are associated with spinelessness. (B) Manhattan plot showing genetic differentiation on chromosome 9 between spined and spineless groups from Lake Storvatnet, using a sliding window genome-wide Fisher’s exact test. A highly differentiated “F_ST-island” located between 4 and 4.4 Mb. (C) Nucleotide diversity (π) and Tajima’s D plots for the 4 to 4.4 Mb region. The sky-blue line represents the spineless group, which shows higher genetic diversity than spined group (red line) in this region.

3.3.2 The Pitx1 locus

Approximately 11.41 M, 10.87 M, 11.17 M, and 11.46 M paired-end reads were successfully mapped for Spined-Storvatnet, Spineless-Storvatnet, Gjerdhaugvatnet, and marine PoolSeq data, respectively, against the Pitx1 BAC library of the threespine stickleback (GenBank: GU130435.1), with an average mapping quality of 35.6 across all population samples. The median sequencing depth across the Pitx1 locus was 225x.

While analyzing the PoolSeq coverage depth along the BAC Pitx1 locus, we identified a notable reduction in depth between positions 129973 to 130025 (~52 bp), located adjacent to the TG-III repeats within the PelA region, in samples from Lake Storvatnet (both spined and spineless) compared to those from Lake Gjerdhaugvatnet and the marine site. In the Storvatnet specimens, this reduction in coverage depth indicates the presence of two alleles, with one allele carrying the deletion (~52 bp) and the other without. Interestingly, the coverage depth in the PoolSeq data for the spineless group is not zero, but remains lower than in the spined group, suggesting that a higher proportion of individuals in the spineless group carry the 52 bp deletion in both alleles. This pattern points to the deletion being more prevalent among the spineless specimens (Figures 9A–C).

Figure 9

Figure 9. Analysis of the PelA region across stickleback populations using data from 40 pooled individuals and from nanopore sequencing. (A) Coverage depth plots for the PelA region (BAC sequence) using pooled sequencing, comparing populations from Lake Gjerdhaugvatnet (blue), the marine environment (gray), spined individuals from Lake Storvatnet (red), and spineless individuals from Lake Storvatnet (light blue). Differences in sequencing coverage indicate structural variation, with a deletion observed in the spineless group near position 130,000 bp. (B) Nanopore sequencing depth for the PelA region from one specimen, highlighting the absence of aligned sequences near position 130,000 bp in the spineless group, supporting that a deletion has occurred in this region. One spined and one spineless specimen from Lake Storvatnet were subject to Nanopore sequencing. (C) DNA sequence of the PelA enhancer region, showing the location of TG repeats. Red sequences inside highlighted area indicate the deleted sequences (52 bp) adjacent to the TG-III repeats in spineless groups (D) F_ST plot across the PelA region between spined and spineless groups from Lake Storvatnet, showing a differentiation peak around position 130,000 bp, where deletion is suggested. (E) Nucleotide diversity (π) and Tajima’s D plots showing differences between spined and spineless individuals in the PelA region.

Additionally, individual nanopore duplex sequencing of a spined and a spineless stickleback from Lake Storvatnet produced 191,722, and 144,670 mapped reads, respectively, with an average mapping quality of 40 for each individual. The median nanopore-sequencing-coverage-depths for the spined and spineless specimens were 18x and 11x, respectively. Nanopore sequencing of the spineless individual showed zero coverage depth in the 52 bp region, consistent with a deletion at both alleles, while the spined individual exhibited higher coverage depth in this region, consistent with the pooled sequencing data (Figures 9A–C). These findings suggest that the 52 bp deletion is present in the spineless individual, most likely in the homozygous state, whereas in the spined individual, the deletion is either absent or it occurs in a heterozygous state.

When assessing -log₁₀(p) values between spined and spineless groups along this region, we observed elevated peaks, suggesting some degree of genetic differentiation. Additionally, the analysis of nucleotide diversity and T_D revealed that the spined group has slightly higher nucleotide diversity and more negative but higher T_D values, which is indicative of balancing selection (Figures 9D, E). These data suggest that the spined group contains more heterozygous individuals carrying both the deleted and non-deleted alleles, while the spineless group has a higher proportion of homozygous individuals with the deletion in both alleles.

4 Discussion

Genomic analysis of marine and freshwater sticklebacks revealed high genetic diversity among marine sticklebacks while two freshwater populations that are located close to each other in the same watercourse are genetically more different from each other than they are from their marine conspecifics. We identified genomic regions with signatures of selection, contributing to our understanding of evolutionary forces shaping these populations. We also identified variation in Pitx1 enhancer regions and other genomics regions that might contribute to pelvic spine reduction.

4.1 Technical implications of pooled DNA sequencing

Although the PoolSeq method has limitations, such as the loss of information on haplotypes, heterozygosity, and linkage disequilibrium (Cutler and Jensen, 2010), it remains a cost-effective strategy for collecting representative data on population samples (Anand et al., 2016). Typically, DNA pools with ≥30 individual DNA samples offer a reliable estimate of allele frequencies (Gautier et al., 2013; Rode et al., 2018). We used a sample size of 40 individuals per DNA pool coupled with high-coverage and deep sequencing to improve allele frequency estimation and decrease the risk of false positives (Cutler and Jensen, 2010; Ferretti et al., 2013; Gautier et al., 2013; Rode et al., 2018). To further validate our findings, we supplemented the PoolSeq data with nanopore long-read sequencing of one spined and one spineless stickleback from Lake Storvatnet, both of which were included in the PoolSeq samples. Nanopore sequencing, known for its ability to read long DNA fragments, typically exceeding 10 kb, is particularly effective at resolving repetitive DNA sequences. This approach enabled a more accurate reconstruction of complex genomic regions, providing clearer and more comprehensive insights into the DNA sequence (Amarasinghe et al., 2020).

4.2 Genetic diversity and differentiation between populations

The upper Lake Storvatnet (altitude 80 m) and Lake Gjerdhaugvatnet (altitude 20 m) exhibit different levels of genetic diversity compared to marine sticklebacks, consistent with Klepaker et al. (2012). Recolonization of Northern Europe by threespine sticklebacks began approximately 17.1 – 37.3 thousand years ago (Fang et al., 2018). As the freshwater lakes like Lake Gjerdhaugvatnet and Lake Storvatnet were likely colonized by small subsets of the marine stickleback populations, this have resulted in a founder effect, leading to reduced genetic diversity in these lakes. Genetic diversity is influenced by the interplay of mutations, gene flow, genetic drift and natural selection. While mutations are rare and selection targets specific genetic loci, overall diversity is largely determined by the long-term effective population size (N_e), which is conversely proportional to genetic drift. Marine sticklebacks, with their moderate dispersal capabilities, possibly reinforced by site fidelity, likely exhibit some degree of population structuring. However, their large census population sizes (N) and the lack of impassable barriers to gene flow in the marine environment likely result in high effective population sizes, translating into less genetic drift and higher overall genetic diversity compared to freshwater populations. Since only a single marine site, a river outlet, was sampled in this study, the full extent of marine genetic diversity remains uncertain.

Sticklebacks in Lake Storvatnet and Lake Gjerdhaugvatnet are highly differentiated, a finding supported by the signatures of selection observed in specific genomic regions of several autosomes. The presence of both spined and pelvic reduced sticklebacks in Lake Storvatnet further underscores this differentiation. Firstly, the colonization history of these lakes is unclear, but given that “the marine limit” (the maximum altitude of the sea surface relative to today’s sea level since the last ice-age; Geological Survey of Norway) in this area is approximately 35 m above sea level, [Available online at: https://geo.ngu.no/kart/losmasse_mobil/?lang=nor&map=9 (Accessed February 20, 2025)] sticklebacks could have reached Lake Storvatnet (80 m altitude) only through non-typical means, such as transport by humans or birds, natural disasters such as tsunami, or via ancient streams that are no longer present. In contrast, Lake Gjerdhaugvatnet, at 20 m, sits below the marine limit, and may have experienced more frequent and prolonged gene flow from nearby marine stickleback populations, influencing its genetic composition. Secondly, genetic drift may have a greater impact on the population in Lake Gjerdhaugvatnet due to its smaller size and potential for population fluctuations caused by environmental factors like weather. Although population size estimates are unavailable, it is reasonable to assume that Lake Storvatnet, being 8–10 times larger, supports a higher effective population size (Ne), which would make it less susceptible to drift. In contrast, the smaller Lake Gjerdhaugvatnet may experience more pronounced effects of drift, leading to reduced genetic diversity. Thirdly, if these differences in population size and drift hold, it could render natural selection on adaptive traits, such as pelvic spine morphology, more effective in Lake Storvatnet. The genomic evidence, particularly the distinct signals of directional selection in Lake Storvatnet (e.g., regions on Chr 1, 7, 9, and 17), underscores the role of selection as an important evolutionary force shaping this subarctic freshwater stickleback population.

4.3 Genes potentially involved in bone or limb development

Based on population genetic differentiation among sites we identified several genes that might contribute to the adaptive structuring of stickleback populations. Genes such as Tbx4 and Fgf8a are known to be crucial in pelvic fin and limb development, with their reduced expression linked to pelvic reduction in some fish species (Cole et al., 2003; Petit et al., 2017). Similarly, Pou1f1 and Alx homeobox genes play essential roles in growth regulation, pituitary organogenesis, and limb and craniofacial development (Cole et al., 2003; Işık and Bilgen, 2019; McGonnell et al., 2011; Petit et al., 2017). Astn2, Brinp, Eng, and Lhx3 are involved in embryonic development, neuronal identity, and brain development in vertebrates (Berkowicz et al., 2016; Srivastava et al., 2010). Genes Gdf10, Bmp2, and Sox9a play vital roles in interdigital webbing, limb morphogenesis, dorsoventral patterning, inhibition of osteoblast differentiation, and cartilage formation (Cheng et al., 2016; Cresko et al., 2003; Kishimoto et al., 1997; Lai et al., 2006). Similarly, genes like H2az, Mllt3, Bbx, Bsx, and H2a have broad implications in chromatin structure, hematopoietic stem cell maintenance, central nervous system development, and brain-specific functions (Chen et al., 2014; Cremona et al., 2004; Germano et al., 2022; Giaimo et al., 2019).

The above-mentioned genes were identified mainly based on population genetic differentiation of the genomic region in which they reside, which is influenced by demographic factors such as founder effects and genetic drift as well. Also, the present study lacks the resolution, including data on recombination frequencies, needed to conclude on the effect of specific loci. Thus, the evidence of candidate genes remains inconclusive, nevertheless, several genes such as H2az and Mllt3 on chromosome 7 (see Supplementary S14), are located within genomic regions displaying patterns in line with selection, potentially reflecting selective pressures acting on these regions and linked to functions such as limb or bone development.

4.4 Parallel evolution of chromosome 4 in freshwaters

Chromosome 4 shows notable genetic differentiation between freshwater and marine stickleback populations, underscoring its importance in adaptation to freshwater environments. Several regions on this chromosome, including the Eda locus, exhibit pronounced genetic differentiation, suggesting strong selection pressures driving divergence between marine and freshwater populations. In our study, the absence of coverage-depth in intronic regions near exon 2 of the Eda locus in freshwater populations, along with SNPs showing high allele frequency differences between freshwater and marine populations, suggests an interplay between structural variation and allelic diversity contributing to the genetic divergence between these populations. The presence of fixed SNPs in intergenic regions, including a fixed G base at position 12,794,969 (Table 2) between the Eda and Tnfsf13b loci in freshwater sticklebacks and a polymorphic (G/T) base at the same position in marine sticklebacks, aligns with previous findings (Jones et al., 2012; Laurentino et al., 2022; O’Brown et al., 2015; Rodríguez-Ramírez et al., 2023) that differentiate low-plated from fully plated sticklebacks. These patterns support the concept of parallel local adaptation, where similar selective forces have repeatedly shaped genetic variation in freshwater sticklebacks, potentially influenced by standing genetic variation (SGV) present in this chromosomal region. Moreover, the absence of coverage-depth in intronic regions near exon 2 of the Eda locus in freshwater populations may indicate structural modifications at this locus, which could play a role in shaping the evolutionary trajectory of chromosome 4 and contribute to the distinct adaptations observed in these subarctic freshwater populations.

4.5 Role of chromosomal inversion in local adaptation

We detected three major putative chromosomal inversions of sizes 0.5 Mb, 0.4 Mb, and 2 Mb on chromosomes 1, 11, and 21, respectively. These inversions are likely to be associated with freshwater adaptations. Similar genetic adaptations have been reported in stickleback populations from North- America and Russia (Jones et al., 2012; Terekhanova et al., 2014). Jones et al. (2012) also provide evidence for the role of chromosomal inversions in divergent selection of sticklebacks in marine and freshwater environments. Hence, these repeated genome changes seem to be examples of the reuse of shared SGV at a global scale (Jones et al., 2012).

In recent years inversions have gained attention for their role in protecting large chromosomal regions containing hundreds of genes from meiotic recombination, thus facilitating local adaptation. Genomic studies have demonstrated that inversions are taxonomically diverse, and often large in size, ranging from 130 Kb to 100 Mb across plants and animals (reviewed by Wellenreuther and Bernatchez, 2018). Inversions are often old and maintained through balancing selection, which indicate their role in maintaining genetic diversity (Wellenreuther, 2017). They also have broad biological impacts, including phenotype–genotype associations, mating behavior, and environmental adaptation. Moreover, the varying proportions of inversions between our two freshwater lakes might reflect reproductive isolation and genetic drift, as previously suggested (Klepaker et al., 2012).

4.6 Divergence between spined and spineless sticklebacks within Lake Storvatnet

4.6.1 Chromosome 9

The Hand2 gene, located within the “F_ST-island” on chr 9, is known to be involved in pathways related to pelvis development in vertebrates (reviewed in Swank et al., 2021). The Hand2 also plays a role in diverse cellular processes related to heart and forelimb development in zebrafish and mouse (Osterwalder et al., 2014; Yelon et al., 2000). Along with its upstream enhancer, Hand2 collaborates with Hoxd13 to activate the Sonic Hedgehog (SHH) pathway during limb bud development (Galli et al., 2010). This SHH pathway is crucial for limb patterning, as it inhibits the GLI3 repressor (GLI3R), preventing it from suppressing Hand2 activity (see Supplementary S9) (Osterwalder et al., 2014; Zuniga, 2015). Disruptions in this pathway can lead to a loss of limb polarity and defects in limb development (Galli et al., 2010).

In addition, the spineless group exhibits higher π and elevated T_D values compared to the spined group within the “F_ST-island” on chromosome 9, suggesting that balancing selection may be maintaining heterozygous alleles at higher frequencies in the spineless group. The persistence of heterozygous alleles in these regions could contribute to the spineless phenotype by enabling the population to maintain genetic variation essential for this trait. Given these insights, our study opens up avenues for further research into the role of Hand2 in the development of pelvic spines in threespine sticklebacks, especially since there are currently no reported studies demonstrating the involvement of Hand2 in pelvic spine reduction.

4.6.2 Pitx1

The Pitx1 locus, situated at the sub-telomeric end of chromosome 7, contains repetitive DNA characterized by TG dinucleotide repeats, which makes sequence assembly difficult (Chan et al., 2010). This inherent complexity likely contributes to its absence from the latest reference genome assemblies. Previous studies (Chan et al., 2010; Thompson et al., 2018) have linked indels at the PelA and PelB enhancers of the Pitx1 locus to reduced pelvic spine length in North American sticklebacks. In our prior work, we identified deletions adjacent to TG-III repeats within the PelA enhancer in both spined and spineless sticklebacks from Lake Storvatnet, which were absent in marine and Lake Gjerdhaugvatnet populations, but found no clear association with pelvic spine reduction (Adhikari et al., 2023).

Expanding on this, our current study suggests that deletions near the TG-III repeats are more prevalent in spineless sticklebacks from Lake Storvatnet, where they predominantly occur as homozygous deletions. This homozygous state means that both copies of the gene have the deletion, likely acting as a recessive allele favoring the spineless phenotype. In contrast, the spined group appears to exhibit a higher frequency of the heterozygous genotype, where one copy of the gene has the deletion and the other does not. This heterozygote state suggests a dominance of non-deleted allele, possibly favoring pelvic spine development. This pattern was consistently observed in both PoolSeq and nanopore sequence data analyses. Further, F_ST values indicated some degree of genetic differentiation between spined and spineless groups in this region. The slightly higher nucleotide diversity and more negative, but higher T_D values in the spined group compared to the spineless group, suggest that it likely contains more heterozygous individuals, indicating balancing selection. In contrast, the spineless group showed a higher proportion of homozygous individuals, potentially linking this genetic pattern to the spineless phenotype by influencing pelvic spine development.

These observations are in line with the concept that recessive alleles, such as those leading to pelvic reduction, may be carried at low frequencies in the spined individuals (heterozygous state), while the spineless phenotype is more likely to be expressed in individuals that are homozygous for the deletion. Thus, the presence of both homozygous and heterozygous states in these populations underscores the role of recessive and dominant genetic mechanisms in driving phenotypic evolution in sticklebacks (Bell et al., 2007).

As the PelA deletion observed in our study was smaller than those documented in North American sticklebacks, we speculate that these deletions might be at a tipping point for transcriptional activity in Lake Storvatnet sticklebacks. It is therefore likely that loci other than PelA and PelB, such as Hand2 and others, contribute to the observed pelvic spine reduction, pointing to a more complex genetic architecture than previously thought. Such findings suggest a case of parallel evolution, where similar phenotypic traits arise independently through different genetic mutations. A similar phenomenon has been observed in Peromyscus polionotus populations in Florida, where variations in fur color arose through different genetic changes (Hoekstra et al., 2006). Indeed, as noted by Poore et al. (2022), it remains uncertain how frequently genetic parallelism, i.e. the use of the same genetic changes, underlies phenotypic parallelism in sticklebacks.

5 Conclusion

In this study, we present evidence of genomic differentiation between marine and freshwater sticklebacks in the European subarctic, and the genomic signatures of local adaptation in two closely located freshwater stickleback populations with differing pelvic morphologies. The study revealed genetic differences between spined and spineless sticklebacks within the Lake Storvatnet. Our previous study using Sanger sequencing found no clear association between deletions in Pitx1 enhancers and pelvic reduction in Lake Storvatnet. The present study, based on deep sequencing of pooled population samples and long-read sequencing of two specimens suggests that homozygous deletions adjacent to TG-III repeats within the PelA region may contribute to this trait. Validation with a larger sample of individual Nanopore long-read sequences of spined and spineless individuals could confirm or disprove these findings. Additionally, we identified a differentiated “F_ST-island” on chromosome 9 between the spined and spineless groups, which includes Hand2, a gene known to be involved in limb development and potentially playing a supportive role in pelvic spine reduction in the Lake Storvatnet population. Together, these findings indicate that pelvic reduction in these sticklebacks is caused by the interplay of several genes, suggesting a case of parallel evolution shaped by polygenic effects rather than a single-gene effect.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, Bioproject PRJNA1073316.

Ethics statement

The animal study was approved by Norwegian Ministry of Agriculture and Food. The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

DA: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. BK: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – review & editing. TJ: Methodology, Writing – review & editing, Investigation. SJ: Conceptualization, Investigation, Methodology, Supervision, Validation, Writing – review & editing. JN: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Writing – review & editing. TM: Conceptualization, Data curation, Investigation, Methodology, Supervision, Validation, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was supported by general grants from Nord University.

Acknowledgments

We thank the genomic facility at Nord University for general support. We also thank Karoline Lauritzen, a former MSc student at Nord University, for helping with the Nanopore sequencing.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2025.1546874/full#supplementary-material

References

Available online at: https://geo.ngu.no/kart/losmasse_mobil/?lang=nor&map=9 (Accessed February 20, 2025).

Google Scholar

Adhikari D., Hanssen I. K., Johansen S. D., Moum T. B., Nordeide J. T. (2023). Pitx1 enhancer variants in spined and spine-reduced subarctic European sticklebacks. Fishes 8, 164. doi: 10.3390/fishes8030164

Crossref Full Text | Google Scholar

Amarasinghe S. L., Su S., Dong X., Zappia L., Ritchie M. E., Gouil Q. (2020). Opportunities and challenges in long-read sequencing data analysis. Genome Biol. 21, 30. doi: 10.1186/s13059-020-1935-5

PubMed Abstract | Crossref Full Text | Google Scholar

Anand S., Mangano E., Barizzone N., Bordoni R., Sorosina M., Clarelli F., et al. (2016). Next Generation Sequencing of pooled samples: Guideline for variants’ filtering. Sci. Rep. 6, 33735. doi: 10.1038/srep33735

PubMed Abstract | Crossref Full Text | Google Scholar

Arendt J., Reznick D. (2008). Convergence and parallelism reconsidered: what have we learned about the genetics of adaptation? Trends Ecol. Evol. 23, 26–32. doi: 10.1016/j.tree.2007.09.011

PubMed Abstract | Crossref Full Text | Google Scholar

Barrett R. D. H., Schluter D. (2008). Adaptation from standing genetic variation. Trends Ecol. Evol. 23, 38–44. doi: 10.1016/j.tree.2007.09.008

PubMed Abstract | Crossref Full Text | Google Scholar

Bell M. A. (1987). Interacting evolutionary constraints in pelvic reduction of threespine sticklebacks, Gasterosteus aculeatus (Pisces, Gasterosteidae). Biol. J. Linn Soc 31, 347–382. doi: 10.1111/j.1095-8312.1987.tb01998.x

Crossref Full Text | Google Scholar

Bell M. A. (1994). Palaeobiology and evolution of the threespine stickleback” in The evolutionary biology of the threespine stickleback. Eds. Bell M. A., Foster S. A. (Oxford: Oxford University Press), 438–471.

Google Scholar

Bell M. A., Khalef V., Travis M. P. (2007). Directional asymmetry of pelvic vestiges in threespine stickleback. J. Exp. Zool. B: Mol. Dev. Evol. 308B, 189–199. doi: 10.1002/jez.b.21132

PubMed Abstract | Crossref Full Text | Google Scholar

Bell M. A., Ortí G. (1994). Pelvic reduction in threespine stickleback from Cook Inlet lakes: geographical distribution and intrapopulation variation. Copeia 1994(2), 314–325. doi: 10.2307/1446981

Crossref Full Text | Google Scholar

Bell M. A., Orti G., Walker J. A., Koenings J. P. (1993). Evolution of pelvic reduction in threespine stickleback fish: a test of competing hypotheses. Evolution 47, 906–914. doi: 10.1111/j.1558-5646.1993.tb01243.x

PubMed Abstract | Crossref Full Text | Google Scholar

Berkowicz S. R., Featherby T. J., Qu Z., Giousoh A., Borg N. A., Heng J. I., et al. (2016). Brinp1^–/– mice exhibit autism-like behaviour, altered memory, hyperactivity and increased parvalbumin-positive cortical interneuron density. Mol. Autism. 7, 22. doi: 10.1186/s13229-016-0079-7

PubMed Abstract | Crossref Full Text | Google Scholar

Chan Y. F., Marks M. E., Jones F. C., Villarreal G., Shapiro M. D., Brady S. D., et al. (2010). Adaptive evolution of pelvic reduction in sticklebacks by recurrent deletion of a Pitx1 Enhancer. Science 327, 302–305. doi: 10.1126/science.1182213

PubMed Abstract | Crossref Full Text | Google Scholar

Chen T., Zhou L., Yuan Y., Fang Y., Guo Y., Huang H., et al. (2014). Characterization of Bbx, a member of a novel subfamily of the HMG-box superfamily together with Cic. Dev. Genes Evol. 224, 261–268. doi: 10.1007/s00427-014-0476-x

PubMed Abstract | Crossref Full Text | Google Scholar

Cheng C. W., Hsiao J. R., Fan C. C., Lo Y. K., Tzen C. Y., Wu L. W., et al. (2016). Loss of GDF10/BMP3b as a prognostic marker collaborates with TGFBR3 to enhance chemotherapy resistance and epithelial-mesenchymal transition in oral squamous cell carcinoma. Mol. Carcinog. 55, 499–513. doi: 10.1002/mc.22297

PubMed Abstract | Crossref Full Text | Google Scholar

Cole N. J., Tanaka M., Prescott A., Tickle C. (2003). Expression of limb initiation genes and clues to the morphological diversification of threespine stickleback. Curr. Biol. 13, R951–R952. doi: 10.1016/j.cub.2003.11.039

PubMed Abstract | Crossref Full Text | Google Scholar

Colosimo P. F., Hosemann K. E., Balabhadra S., Villarreal G., Dickson M., Grimwood J., et al. (2005). Widespread parallel evolution in sticklebacks by repeated fixation of Ectodysplasin alleles. Science 307, 1928–1933. doi: 10.1126/science.1107239

PubMed Abstract | Crossref Full Text | Google Scholar

Colosimo P. F., Peichel C. L., Nereng K., Blackman B. K., Shapiro M. D., Schluter D., et al. (2004). The genetic architecture of parallel armor plate reduction in threespine sticklebacks. PLoS Biol. 2, e109. doi: 10.1371/journal.pbio.0020109

PubMed Abstract | Crossref Full Text | Google Scholar

Coyle S. M., Huntingford F. A., Peichel C. L. (2007). Parallel evolution of Pitx1 underlies pelvic reduction in Scottish threespine stickleback (Gasterosteus aculeatus). J. Hered. 98, 581–586. doi: 10.1093/jhered/esm066

PubMed Abstract | Crossref Full Text | Google Scholar

Cremona M., Colombo E., Andreazzoli M., Cossu G., Broccoli V. (2004). Bsx, an evolutionary conserved Brain Specific homeoboX gene expressed in the septum, epiphysis, mammillary bodies and arcuate nucleus. Gene Expr. Patterns. 4, 47–51. doi: 10.1016/S1567-133X(03)00151-0

PubMed Abstract | Crossref Full Text | Google Scholar

Cresko W. A., Amores A., Wilson C., Murphy J., Currey M., Phillips P., et al. (2004). Parallel genetic basis for repeated evolution of armor loss in Alaskan threespine stickleback populations. Proc. Natl. Acad. Sci. U.S.A. 101, 6050–6055. doi: 10.1073/pnas.0308479101

PubMed Abstract | Crossref Full Text | Google Scholar

Cresko W. A., Yan Y.-L., Baltrus D. A., Amores A., Singer A., Rodríguez-Marí A., et al. (2003). Genome duplication, subfunction partitioning, and lineage divergence: Sox9 in stickleback and zebrafish. Dev. Dynam. 228, 480–489. doi: 10.1002/dvdy.10424

PubMed Abstract | Crossref Full Text | Google Scholar

Cutler D. J., Jensen J. D. (2010). To pool, or not to pool? Genetics 186, 41–43. doi: 10.1534/genetics.110.121012

PubMed Abstract | Crossref Full Text | Google Scholar

DeFaveri J., Shikano T., Shimada Y., Goto A., Merilä J. (2011). Global analysis of genes involved of genes involved in freshwater adaptation in threespine sticklebacks (Gasterosteus aculeatus). Evolution 65, 1800–1807. doi: 10.1111/j.1558-5646.2011.01247.x

PubMed Abstract | Crossref Full Text | Google Scholar

Fang B., Kemppainen P., Momigliano P., Feng X., Merilä J. (2020). On the causes of geographically heterogeneous parallel evolution in sticklebacks. Nat. Ecol. Evol. 4, 1105–1115. doi: 10.1038/s41559-020-1222-6

PubMed Abstract | Crossref Full Text | Google Scholar

Fang B., Merilä J., Ribeiro F., Alexandre C. M., Momigliano P. (2018). Worldwide phylogeny of three-spined sticklebacks. Mol. Phylogenet. Evol. 127, 613–625. doi: 10.1016/j.ympev.2018.06.008

PubMed Abstract | Crossref Full Text | Google Scholar

Ferretti L., Ramos-Onsins S. E., Pérez-Enciso M. (2013). Population genomics from pool sequencing. Mol. Ecol. 22, 5561–5576. doi: 10.1111/mec.12522

PubMed Abstract | Crossref Full Text | Google Scholar

Galli A., Robay D., Osterwalder M., Bao X., Bénazet J.-D., Tariq M., et al. (2010). Distinct roles of Hand2 in initiating polarity and posterior Shh expression during the onset of mouse limb bud development. PLoS Genet. 6, e1000901. doi: 10.1371/journal.pgen.1000901

PubMed Abstract | Crossref Full Text | Google Scholar

Gautier M., Foucaud J., Gharbi K., Cézard T., Galan M., Loiseau A., et al. (2013). Estimation of population allele frequencies from next-generation sequencing data: pool-versus individual-based genotyping. Mol. Ecol. 22, 3766–3779. doi: 10.1111/mec.12360

PubMed Abstract | Crossref Full Text | Google Scholar

Germano G., Porazzi P., Felix C. A. (2022). Leukemia-associated transcription factor mllt3 is important for primitive erythroid development in zebrafish embryogenesis. Dev. Dynam. 251, 1728–1740. doi: 10.1002/dvdy.477

PubMed Abstract | Crossref Full Text | Google Scholar

Giaimo B. D., Ferrante F., Herchenröther A., Hake S. B., Borggrefe T. (2019). The histone variant H2A.Z in gene regulation. Epigenet. Chromatin. 12, 37. doi: 10.1186/s13072-019-0274-9

PubMed Abstract | Crossref Full Text | Google Scholar

Gregorius H. R. (1987). The relationship between the concepts of genetic diversity and differentiation. Theor. Appl. Genet. 74, 397–401. doi: 10.1007/BF00274724

PubMed Abstract | Crossref Full Text | Google Scholar

Hagen D. W., Gilbertson L. G. (1973). Selective predation and the intensity of selection acting upon the lateral plates of threespine sticklebacks. Heredity 30, 273–287. doi: 10.1038/hdy.1973.38

Crossref Full Text | Google Scholar

Hoekstra H. E., Hirschmann R. J., Bundey R. A., Insel P. A., Crossland J. P. (2006). A single amino acid mutation contributes to adaptive beach mouse color pattern. Science 313, 101–104. doi: 10.1126/science.1126121

PubMed Abstract | Crossref Full Text | Google Scholar

Hohenlohe P. A., Bassham S., Etter P. D., Stiffler N., Johnson E. A., Cresko W. A. (2010). Population genomics of parallel adaptation in threespine stickleback using sequenced RAD tags. PloS Genet. 6, e1000862. doi: 10.1371/journal.pgen.1000862

PubMed Abstract | Crossref Full Text | Google Scholar

Innan H., Kim Y. (2004). Pattern of polymorphism after strong artificial selection in a domestication event. Proc. Natl. Acad. Sci. U.S.A. 101, 10667–10672. doi: 10.1073/pnas.0401720101

PubMed Abstract | Crossref Full Text | Google Scholar

Işık R., Bilgen G. (2019). Associations between genetic variants of the POU1F1 gene and production traits in Saanen goats. Arch. Anim. Breed. 62, 249–255. doi: 10.5194/aab-62-249-2019

PubMed Abstract | Crossref Full Text | Google Scholar

Jones F. C., Grabherr M. G., Chan Y. F., Russell P., Mauceli E., Johnson J., et al. (2012). The genomic basis of adaptive evolution in threespine sticklebacks. Nature 484, 55–61. doi: 10.1038/nature10944

PubMed Abstract | Crossref Full Text | Google Scholar

Karlsen B. O., Adhikari D., Jørgensen T. E., Hanssen I. K., Moum T. B., Nordeide J. T., et al. (2024). Two distinct maternal lineages of threespine stickleback (Gasterosteus aculeatus) in a small Norwegian subarctic lake. Fishes 9, 285. doi: 10.3390/fishes9070285

Crossref Full Text | Google Scholar

Kishimoto Y., Lee K.-H., Zon L., Hammerschmidt M., Schulte-Merker S. (1997). The molecular nature of zebrafish swirl: BMP2 function is essential during early dorsoventral patterning. Development 124, 4457–4466. doi: 10.1242/dev.124.22.4457

PubMed Abstract | Crossref Full Text | Google Scholar

Klepaker T. (1993). Morphological-changes in a marine population of threespine tickleback, Gasterosteus aculeatus, recently isolated in freshwater. Can. J. Zool. 71, 1251–1258. doi: 10.1139/z93-171

Crossref Full Text | Google Scholar

Klepaker T. (1995). Postglacial evolution in lateral plate morphs in Norwegian freshwater populations of the threespine stickleback (Gasterosteus aculeatus). Can. J. Zool. 73, 898–906. doi: 10.1139/z95-105

Crossref Full Text | Google Scholar

Klepaker T. O., Østbye K. (2008). Pelvic anti-predator armour reduction in Norwegian populations of the threespine stickleback: a rare phenomenon with adaptive implications? J. Zool. 276, 81–88. doi: 10.1111/j.1469-7998.2008.00471.x

Crossref Full Text | Google Scholar

Klepaker T., Ostbye K., Bell M. A. (2013). Regressive evolution of the pelvic complex in stickleback fishes: a study of convergent evolution. Evol. Ecol. Res. 15, 413–435.

Google Scholar

Klepaker T., Ostbye K., Bernatchez L., Vøllestad L. A. (2012). Spatio-temporal patterns in pelvic reduction in threespine stickleback (Gasterosteus aculeatus L.) in Lake Storvatnet. Evol. Ecol. Res. 14, 169–191.

Google Scholar

Kofler R., Orozco-terWengel P., De Maio N., Pandey R., Nolte V., Futschik A., et al. (2011a). PoPoolation: A toolbox for population genetic analysis of next generation sequencing data from pooled individuals. PLoS One 6, e15925. doi: 10.1371/journal.pone.0015925

PubMed Abstract | Crossref Full Text | Google Scholar

Kofler R., Pandey R., Schlötterer C. (2011b). PoPoolation2: identifying differentiation between populations using sequencing of pooled DNA samples (Pool-Seq). Bioinformatics 27, 3435–3436. doi: 10.1093/bioinformatics/btr589

PubMed Abstract | Crossref Full Text | Google Scholar

Lai C.-F., Bai S., Uthgenannt B. A., Halstead L. R., McLoughlin P., Schafer B. W., et al. (2006). Four and half lim protein 2 (FHL2) stimulates osteoblast differentiation. J. Bone Miner Res. 21, 17–28. doi: 10.1359/JBMR.050915

PubMed Abstract | Crossref Full Text | Google Scholar

Laurentino T. G., Boileau N., Ronco F., Berner D. (2022). The ectodysplasin-A receptor is a candidate gene for lateral plate number variation in stickleback fish. G3: Genes Genomes Genet. 12(6). doi: 10.1093/g3journal/jkac077

PubMed Abstract | Crossref Full Text | Google Scholar

Liu C., Fetterman J. L., Liu P., Luo Y., Larson M. G., Vasan R. S., et al. (2018). Deep sequencing of the mitochondrial genome reveals common heteroplasmic sites in NADH dehydrogenase genes. Hum. Genet. 137, 203–213. doi: 10.1007/s00439-018-1873-4

PubMed Abstract | Crossref Full Text | Google Scholar

McGonnell I. M., Graham A., Richardson J., Fish J. L., Depew M. J., Dee C. T., et al. (2011). Evolution of the Alx homeobox gene family: parallel retention and independent loss of the vertebrate Alx3 gene. Evol. Dev. 13, 343–351. doi: 10.1111/j.1525-142X.2011.00489.x

PubMed Abstract | Crossref Full Text | Google Scholar

McPhail J. D. (1992). Ecology and evolution of sympatric sticklebacks (Gasterosteus): evidence for a species-pair in Paxton Lake, Texada Island, British Columbia. Can. J. Zool. 70, 361–369. doi: 10.1139/z92-054

Crossref Full Text | Google Scholar

Moodie G. E. E., Reimchen T. E. (1976). Phenetic variation and habitat differences in Gasterosteus populations of the Queen Charlotte Islands. Syst. Biol. 25, 49–61. doi: 10.2307/2412778

Crossref Full Text | Google Scholar

Nei M., Li W. H. (1979). Mathematical model for studying genetic variation in terms of restriction endonucleases. Proc. Natl. Acad. Sci. U.S.A. 76, 5269–5273. doi: 10.1073/pnas.76.10.5269

PubMed Abstract | Crossref Full Text | Google Scholar

Nelson T. C., Cresko W. A. (2018). Ancient genomic variation underlies repeated ecological adaptation in young stickleback populations. Evol. Lett. 2, 9–21. doi: 10.1002/evl3.37

PubMed Abstract | Crossref Full Text | Google Scholar

O’Brown N. M., Summers B. R., Jones F. C., Brady S. D., Kingsley D. M. (2015). A recurrent regulatory change underlying altered expression and Wnt response of the stickleback armor plates gene EDA. eLife 4, e05290. doi: 10.7554/eLife.05290

PubMed Abstract | Crossref Full Text | Google Scholar

Osterwalder M., Speziale D., Shoukry M., Mohan R., Ivanek R., Kohler M., et al. (2014). HAND2 targets define a network of transcriptional regulators that compartmentalize the early limb bud mesenchyme. Dev. Cell. 31, 345–357. doi: 10.1016/j.devcel.2014.09.018

PubMed Abstract | Crossref Full Text | Google Scholar

Peichel C. L., Nereng K. S., Ohgi K. A., Cole B. L. E., Colosimo P. F., Buerkle C. A., et al. (2001). The genetic architecture of divergence between threespine stickleback species. Nature 414, 901–905. doi: 10.1038/414901a

PubMed Abstract | Crossref Full Text | Google Scholar

Petit F., Sears K. E., Ahituv N. (2017). Limb development: a paradigm of gene regulation. Nat. Revi. Genet. 18, 245–258. doi: 10.1038/nrg.2016.167

PubMed Abstract | Crossref Full Text | Google Scholar

Poore H. A., Stuart Y. E., Rennison D. J., Roesti M., Hendry A. P., Bolnick D. I., et al. (2022). Repeated genetic divergence plays a minor role in repeated phenotypic divergence of lake-stream stickleback. Evolution 77, 110–122. doi: 10.1093/evolut/qpac025

PubMed Abstract | Crossref Full Text | Google Scholar

Reimchen T. E. (1983). Structural relationships between spines and lateral plates in threespine stickleback (Gasterosteus aculeatus). Evolution 37, 931–946. doi: 10.1111/j.1558-5646.1983.tb05622.x

PubMed Abstract | Crossref Full Text | Google Scholar

Roberts Kingman G. A., Vyas D. N., Jones F. C., Brady S. D., Chen H. I., Reid K., et al. (2021). Predicting future from past: The genomic basis of recurrent and rapid stickleback evolution. Sci. Adv. 7, eabg5285. doi: 10.1126/sciadv.abg5285

PubMed Abstract | Crossref Full Text | Google Scholar

Rode N. O., Holtz Y., Loridon K., Santoni S., Ronfort J., Gay L. (2018). How to optimize the precision of allele and haplotype frequency estimates using pooled-sequencing data. Mol. Ecol. Resour. 18, 194–203. doi: 10.1111/1755-0998.12723

PubMed Abstract | Crossref Full Text | Google Scholar

Rodríguez-Ramírez C. E., Hiltbrunner M., Saladin V., Walker S., Urrutia A., Peichel C. L. (2023). Molecular mechanisms of Eda-mediated adaptation to freshwater in threespine stickleback. Mol. Ecol. 1–19 doi: 10.1111/mec.16989

PubMed Abstract | Crossref Full Text | Google Scholar

Schluter D., Clifford E. A., Nemethy M., McKinnon J. S. (2004). Parallel evolution and inheritance of quantitative traits. Am. Nat. 163, 809–822. doi: 10.1086/383621

PubMed Abstract | Crossref Full Text | Google Scholar

Schluter D., Conte G. L. (2009). Genetics and ecological speciation. Proc. Natl. Acad. Sci. U.S.A. 106, 9955–9962. doi: 10.1073/pnas.0901264106

PubMed Abstract | Crossref Full Text | Google Scholar

Schröder M., Windhager S., Schaefer K., Ahnelt H. (2023). Adaptability of bony armor elements of the threespine stickleback Gasterosteus aculeatus (Teleostei: Gasterosteidae): Ecological and evolutionary insights from symmetry analyses. Symmetry 15, 811. doi: 10.3390/sym15040811

Crossref Full Text | Google Scholar

Sham P., Bader J. S., Craig I., O’Donovan M., Owen M. (2002). DNA Pooling: a tool for large-scale association studies. Nat. Rev. Genet. 3, 862–871. doi: 10.1038/nrg930

PubMed Abstract | Crossref Full Text | Google Scholar

Shapiro M. D., Marks M. E., Peichel C. L., Blackman B. K., Nereng K. S., Jónsson B., et al. (2004). Genetic and developmental basis of evolutionary pelvic reduction in threespine sticklebacks. Nature 428, 717–723. doi: 10.1038/nature02415

PubMed Abstract | Crossref Full Text | Google Scholar

Srivastava M., Larroux C., Lu D. R., Mohanty K., Chapman J., Degnan B. M., et al. (2010). Early evolution of the LIM homeobox gene family. BMC Biol. 8, 4. doi: 10.1186/1741-7007-8-4

PubMed Abstract | Crossref Full Text | Google Scholar

Swank S., Sanger T. J., Stuart Y. E. (2021). (Non)Parallel developmental mechanisms in vertebrate appendage reduction and loss. Ecol. Evol. 11, 15484–15497. doi: 10.1002/ece3.8226

PubMed Abstract | Crossref Full Text | Google Scholar

Tajima F. (1989). Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123, 585–595. doi: 10.1093/genetics/123.3.585

PubMed Abstract | Crossref Full Text | Google Scholar

Terekhanova N. V., Barmintseva A. E., Kondrashov A. S., Bazykin G. A., Mugue N. S. (2019). Architecture of parallel adaptation in ten lacustrine threespine stickleback populations from the White Sea area. Genome Biol. Evol. 11, 2605–2618. doi: 10.1093/gbe/evz175

PubMed Abstract | Crossref Full Text | Google Scholar

Terekhanova N. V., Logacheva M. D., Penin A. A., Neretina T. V., Barmintseva A. E., Bazykin G. A., et al. (2014). Fast evolution from Precast Bricks: Genomics of young freshwater populations of threespine stickleback Gasterosteus aculeatus. PloS Genet. 10, e1004696. doi: 10.1371/journal.pgen.1004696

PubMed Abstract | Crossref Full Text | Google Scholar

Thompson A. C., Capellini T. D., Guenther C. A., Chan Y. F., Infante C. R., Menke D. B., et al. (2018). A novel enhancer near the Pitx1 gene influences development and evolution of pelvic appendages in vertebrates. eLife 7, e38555. doi: 10.7554/eLife.38555

PubMed Abstract | Crossref Full Text | Google Scholar

Wellenreuther M. (2017). Balancing selection maintains cryptic colour morphs. Mol. Ecol. 26, 6185–6188. doi: 10.1111/mec.14406

PubMed Abstract | Crossref Full Text | Google Scholar

Wellenreuther M., Bernatchez L. (2018). Eco-evolutionary genomics of chromosomal inversions. Trends Ecol. Evol. 33, 427–440. doi: 10.1016/j.tree.2018.04.002

PubMed Abstract | Crossref Full Text | Google Scholar

Wootton R. J. (1976). The biology of the sticklebacks (London: Academic Press).

Google Scholar

Yelon D., Ticho B., Halpern M. E., Ruvinsky I., Ho R. K., Silver L. M., et al. (2000). The bHLH transcription factor Hand2 plays parallel roles in zebrafish heart and pectoral fin development. Development 127, 2573–2582. doi: 10.1242/dev.127.12.2573

PubMed Abstract | Crossref Full Text | Google Scholar

Zuniga A. (2015). Next generation limb development and evolution: old questions, new perspectives. Development 142, 3810–3820. doi: 10.1242/dev.125757

PubMed Abstract | Crossref Full Text | Google Scholar