https://orcid.org/0000-0002-7731-8944
Abstract
Recent studies have shown that the repeated evolution of similar phenotypes in response to similar ecological conditions (here “parallel evolution”) often occurs through mutations in the same genes. However, many previous studies have focused on known candidate genes in a limited number of systems. Thus, the question of how often parallel phenotypic evolution is due to parallel genetic changes remains open. Here, we used quantitative trait locus (QTL) mapping in F2 intercrosses between lake and stream threespine stickleback (Gasterosteus aculeatus) from four independent watersheds on Vancouver Island, Canada to determine whether the same QTL underlie divergence in the same phenotypes across, between, and within watersheds. We find few parallel QTL, even in independent crosses from the same watershed or for phenotypes that have diverged in parallel. These findings suggest that different mutations can lead to similar phenotypes. The low genetic repeatability observed in these lake-stream systems contrasts with the higher genetic repeatability observed in other stickleback systems. We speculate that differences in evolutionary history, gene flow, and/or the strength and direction of selection might explain these differences in genetic parallelism and emphasize that more work is needed to move beyond documenting genetic parallelism to identifying the underlying causes.
The repeated evolution of similar phenotypes in independent populations in response to similar environmental conditions is generally used as evidence that phenotypic evolution is driven by natural selection (Endler, 1986; Harvey and Pagel, 1991; Losos, 2011; Schluter, 2000). Repeated evolution has historically been treated as an all-or-nothing phenomenon (Bolnick et al., 2018). However, it is increasingly recognized that repeated evolution is best viewed as a continuum, such that particular phenotypes, populations, or species might exhibit very little repeatability, while others might exhibit very strong repeatability, with variation between the extremes along the continuum (Bolnick et al., 2018; Losos, 2011; Oke et al., 2017). This variation in repeatability provides fertile ground for determining the extent to which natural or sexual selection, demography, evolutionary histories, or genetic constraints shape the outcomes of evolutionary trajectories.
Semantic aside
Such repeated evolution has been referred to as parallel and/or convergent evolution, with debate about how to best distinguish between these terms (e.g., Arendt and Reznick, 2008; Bolnick et al., 2018; Conte et al., 2012; Losos, 2011; Rosenblum et al., 2014; Scotland, 2011; Wake et al., 2011). Here, we use the term “parallel evolution,” which we define as the repeated evolution of similar phenotypes in lineages that are independently derived from a common ancestor (Conte et al., 2012; Simpson, 1961) because it reflects the evolutionary history of our model system (threespine stickleback fish; see below) and is therefore commonly used in this system. We further use the term “parallelism” to refer to the continuum of parallel evolution, from nonparallel to highly parallel evolution, and everything in between.
Genetics of phenotypic parallelism
One possible determinant of parallel evolution is the underlying genetic basis of the phenotypes. If the same genes are used repeatedly when the same phenotypes evolve, it suggests that there could be constraints on the types of genes and mutations that contribute to phenotypic evolution and/or in the availability of standing genetic variation. Alternatively, there could be genotypic redundancy whereby there are multiple genetic solutions to evolve similar phenotypes (Barghi et al., 2020; Láruson et al., 2020). There is now appreciable evidence that the same genes and even the same alleles are often used when the same phenotype evolves (Conte et al., 2012; Martin and Orgogozo, 2013; Stern, 2013). However, many studies have used a candidate gene approach, which asks whether genes known to contribute to the evolution of a given phenotype in one population or species also contribute to the evolution of that phenotype in another population or species. Such candidate gene approaches are subject to publication bias and thus may overestimate genetic parallelism and the role of genetic constraints in phenotypic evolution (Conte et al., 2012).
To take a more agnostic approach to address the question of whether the same or different genes and mutations contribute to parallel phenotypic evolution and adaptation, researchers have used two genome-wide approaches. The first is a population genomics approach that asks whether the same genomic regions show signatures of divergent selection across pairs of populations adapting to similar environments. This approach has now been widely used (reviewed by Fraser and Whiting [2020]), but it has important caveats, such as heterogeneity in recombination rates across the genome, which can make it more likely that the same genomic regions show divergence (Berner and Roesti, 2017; Burri, 2018; Cruickshank and Hahn, 2014; Ravinet et al., 2017; Roesti et al., 2012; Wolf and Ellegren, 2017). Furthermore, population genomic studies usually do not explicitly consider phenotypes. The second is a forward genetic approach that asks whether the same quantitative trait loci (QTL) underlie phenotypes that have evolved in parallel in independent populations or species (e.g., Bainbridge et al., 2020; Blankers et al., 2019; Conte et al., 2015; Erickson et al., 2016; Ferris et al., 2015; Kemppainen et al., 2021). Forward genetic approaches also have caveats, including limitations in the ability to detect QTL of small effect, and a bias toward detecting QTL in regions of low recombination (Beavis, 1998; Noor et al., 2001; Rockman, 2012; Slate, 2013). However, because forward genetic approaches focus on phenotypes, this approach enables an explicit comparison of the extent to which phenotypic parallelism is driven by genetic parallelism.
Phenotypic parallelism in the threespine stickleback lake-stream system
Here, we use QTL mapping to investigate the genetic basis of phenotypic parallelism in lake-stream pairs of threespine stickleback (Gasterosteus aculeatus). These fish have repeatedly colonized freshwater habitats across the Northern hemisphere since the end of the last glacial maximum, approximately 15,000 years ago. Although the marine ancestors of these freshwater populations are relatively homogeneous in phenotype, freshwater sticklebacks show phenotypic divergence from their marine ancestors as well as among extant freshwater populations (Bell and Foster, 1994). In this study, we focus on sticklebacks from freshwater lakes and their adjoining streams, which show repeated patterns of phenotypic evolution in Alaska, British Columbia (Vancouver Island and Haida Gwaii), Iceland, Ireland, northern Germany, and Switzerland (Aguirre, 2009; Berner et al., 2008, 2009, 2010; Hendry and Taylor, 2004; Kaeuffer et al., 2012; Lavin and McPhail, 1993; Lucek et al., 2013, 2014; Moser et al., 2012; Ravinet et al., 2013; Reimchen et al., 1985; Reusch et al., 2001; Stuart et al., 2017). However, parallel divergence is stronger at the regional than the global level, and there is variation in phenotypic parallelism even within a geographic region (Paccard et al., 2020; Stuart et al., 2017). In a study of 16 lake-stream pairs from independent watersheds on Vancouver Island, variation in the multivariate direction of phenotypic parallelism was correlated with the multivariate direction of ecological parallelism, suggesting a role for natural selection in driving the patterns of phenotypic parallelism (Stuart et al., 2017).
The extent to which phenotypic parallelism is driven by genetic parallelism is unknown in the lake-stream system. Several population genomic studies have been conducted on lake-stream systems from both North America and Europe and have found relatively low genomic parallelism (i.e., few shared regions of with signature of divergent selection), at either a local or global scale (Deagle et al., 2012; Feulner et al., 2015; Marques et al., 2016; Paccard et al., 2020; Rennison et al., 2019; Roesti et al., 2012, 2015). However, the genetic basis of specific phenotypic traits that diverge between lake and stream sticklebacks has not been broadly investigated, although at least some of the phenotypic divergence between lake and stream sticklebacks has a genetic basis (Berner et al., 2011; Hendry et al., 2011; Moser et al., 2012; Oke et al., 2016; Sharpe et al., 2008). Only a single QTL mapping study has been performed on a cross between a stream female and a lake male from allopatric populations in Switzerland (Berner et al., 2014), which does not allow us to address the contribution of genetic parallelism to phenotypic parallelism among lake-stream systems. By contrast, at least 18 studies have focused on the genetic basis of phenotypic divergence between marine and freshwater populations or between freshwater benthic and limnetic populations (reviewed in Peichel and Marques, 2017). These studies have revealed a high degree of parallelism at the QTL level (Conte et al., 2015; Erickson et al., 2016), as well as at the level of the genes and genetic changes that underlie specific traits (Peichel and Marques, 2017; Reid et al., 2021). In several cases, the repeated use of the same genetic changes is due to selection on standing variation in the ancestral marine population (Colosimo et al., 2005; Howes et al., 2017; Indjeian et al., 2016; Kitano et al., 2010; Miller et al., 2007), although repeated mutation at the same locus also plays a role in the repeated evolution of pelvic reduction in independent freshwater populations (Chan et al., 2010; Shapiro et al., 2004; Xie et al., 2020). These genetic studies have led to the perception that parallel phenotypic evolution in sticklebacks is largely driven by parallel genetic changes. However, the extent to which genetic parallelism underlies phenotypic parallelism in sticklebacks is unknown outside of a few well-studied systems.
To address this knowledge gap, we investigated the genetic basis of repeated phenotypic evolution in four Vancouver Island lake-stream pairs from independent watersheds (Figure 1A), which have been extensively studied (e.g., Berner et al., 2008, 2009; Hendry and Taylor, 2004; Kaeuffer et al., 2012; Roesti et al., 2012; Stuart et al., 2017). We conducted QTL mapping of 63 morphological traits on at least two independent F2 lake-stream intercrosses from each watershed to address three questions: (1) do the same QTL underlie lake-stream differentiation in the same traits between and among watersheds?; (2) do the same QTL underlie lake-stream differentiation in the same traits within a watershed?; and (3) is there a correlation between phenotypic and genetic parallelism?
QTL mapping in four independent lake-stream pairs. (A) Map of Vancouver Island showing the locations of the four lake-stream pairs (B = Boot; M = Misty; P = Pye; R = Roberts), which are found in independent watersheds. The inset shows a lake and a stream male from Roberts Lake (photos courtesy of Daniel Berner). (B) Illustration for how genetic parallelism within watersheds, between watersheds, and among watersheds was calculated based on the proportion of parallel linkage groups. The example illustrates the nine QTL identified for total plate count. NA indicates that the proportion of shared linkage groups could not be calculated due to lack of QTL in one of the comparisons. QTL = quantitative trait locus.
Materials and methods
Ethics statement
Fish were collected and transported with permission from the British Columbia Ministry of the Environment and the Canadian Department of Fisheries and Oceans (collection permits NA12-84189, SARA 283; transfer permits VI13-84054, VI13-89403). Animal experiments were approved by the Fred Hutchinson Cancer Research Center Institutional Animal Care and Use Committee (protocol 1575).
Collection and crosses
Threespine stickleback fish were collected from the lakes and streams of four watersheds (Boot, Misty, Pye, Roberts; Figure 1A) on Vancouver Island, British Columbia, Canada in May 2013 as described in Stuart et al. (2017).
For each watershed, several independent crosses were made between a lake female and a stream male via in vitro fertilization in the field. These grandparental fish were euthanized in buffered MS-222 and stored in 95% ethanol for future DNA extraction and analysis. The fertilized eggs were transported to the stickleback facility at the Fred Hutchinson Cancer Research Center, where the F1 fish were raised to reproductive maturity. To create the F2 generation, F1 full siblings were crossed via in vitro fertilization. Females and males from two Boot, two Misty, three Pye, and two Roberts F1 families were used to generate six Boot, seven Misty, seven Pye, and five Roberts F2 families (Table S1). The F1 generation was euthanized in buffered MS-222, both pectoral fins were removed and stored in 95% ethanol, and bodies were immediately stored in 10% buffered formalin. The F2 fish were raised until adulthood at approximately 9 months of age and then were euthanized in buffered MS-222. Both pectoral fins were removed and stored in 95% ethanol, and bodies were immediately preserved in 10% buffered formalin for morphological analyses.
Morphological analyses
Formalin-fixed F2 sticklebacks were stained with alizarin red and stored in 37% isopropanol (Peichel et al., 2001). Each F2 stickleback was photographed on its left lateral side, as well as ventrally, with a ruler placed for scale under standardized lighting conditions using a Nikon D90 camera. Prior to photographing, each stickleback was pinned to aid with landmarking following Stuart et al. (2017), with pins at the base of the skull on the dorsal midline, the anterior and posterior insertions of the dorsal fin and anal fin, the anterior edge of the pelvic girdle along the midline, the anterior tip of the ectocoracoid, and the insertions of the pectoral fin (Figure S1).
In addition to lateral and ventral photographs, each F2 was patted dry and weighed to obtain mass, left and right bony armor plate counts were recorded, and measurements were taken via digital calipers for the right and left pelvic spine lengths, and standard length (from the tip of the maxilla to the end of the caudal peduncle along the midline). Gill rakers on the first branchial arch were counted in situ by making an incision just under the right gill flap and cutting up along the base of the chin and through the mouth allowing for the right gill cover to be peeled back and removed. All gill rakers positive for alizarin red were counted.
Using the lateral photographs, body shape measurements were collected using 18 landmarks placed on the left side of the stickleback with tpsDig2 software (v2.10; Rohlf, 2006) (Figure S1). Coordinates from tpsDig were used to calculate the centroid size of each individual fish via Generalized Procrustes Analysis (GPA) in the Geomorph package in R (Adams et al., 2020). The GPA produced x- and y-coordinates corrected for scale, location, and rotation, and these transformed x- and y-coordinates were used as phenotypic traits for QTL mapping.
Additional linear measurements were taken from the lateral (11 traits; Figure S2A) and ventral (10 traits; Figure S2B) photos using ImageJ. These linear traits, as well as mass and pelvic spine lengths, were log transformed and size corrected, following Vamosi (2002). Because each F2 family was raised together in a single tank, size correction was done separately for each F2 family to account for tank effects. Number of plates and gill rakers, as well as standard length and centroid size, were not log transformed or size corrected prior to QTL analyses.
Genotyping
DNA was extracted in 96 well plates from fin clips of the 18 grandparent, 42 F1 (Misty F1 female 1.1 was not included), and 1711 F2 fish (Table S1) using the Promega Wizard SV Genomic DNA Purification kits. Fish were genotyped by double-digest restriction associated DNA sequencing (ddRADseq; Peterson et al., 2012), following Stuart et al. (2017). Each library contained 48 fish, necessitating a total of 37 libraries (two grandparent DNAs and 1 F1 DNA were added to more than one library and DNA from two unrelated fish were also added to one library). Twenty libraries were sequenced across two “rapid read” lanes and eight “high output” lanes, and the remaining 17 libraries were sequenced across eight “high output” lanes using 100 bp, paired end sequencing on an Illumina HiSeq 2500 at the Fred Hutchinson Cancer Research Center Genomics Shared Resource.
Sequence variants were identified using a reference-based bioinformatics pipeline. Demultiplexing was done using Stacks version 1.13 (Catchen et al., 2013). For direct comparison of these results to previous studies of the same populations (Stuart et al., 2017; Rennison et al., 2019), reads were aligned to the original stickleback reference genome version 078 (Jones et al., 2012) using BWA v0.7.7 (Li and Durbin, 2009), and subsequent realignment was done with STAMPY v1.0.23 (Lunter and Goodson, 2011). The GATK v3.3.0 (McKenna et al., 2010) best practices workflow (DePristo et al., 2011) was followed except that the MarkDuplicates step was omitted. RealignTargetCreator and IndelRealigner were used to realign reads around indels, and HaplotypeCaller identified single nucleotide polymorphisms (SNPs) in individuals. Joint genotyping was done across all individuals in each F2 cross from a single watershed using GenotypeGVCFs. The results were written to a VCF file containing all variable sites. This file was filtered for a minimum quality score (20) and depth of coverage (minimum of 8 reads and maximum of 100,000) before downstream analyses.
Linkage map construction and QTL analyses
The sequencing quality of many of the grandparental fish was very poor for unknown reasons. Thus, it was not possible to directly identify SNPs that showed fixed differences between the lake and stream grandparents in any watershed. SNPs that are fixed between the grandparents should also be heterozygous in all genotyped F1s within a watershed, which would enable using the same SNPs for linkage mapping and QTL analyses across all F1 families within a watershed. However, there were very few SNPs that met this criteria in any watershed, suggesting that there are very few fixed differences between lake and stream individuals within a watershed, consistent with previous studies (Rennison et al., 2019; Roesti et al., 2012; Stuart et al., 2017). Thus, SNPs that were heterozygous in all genotyped F1 parents within each F1 family were identified using custom R scripts (Table S2).
All linkage mapping and QTL analyses were performed separately in each F1 family in each of the four watersheds, for a total of nine F1 families (Table S1, Table S2). First, SNP loci with less than 25% missing data across all F2s within an F1 family were extracted to create each linkage map in JoinMap version 4.1 (Van Ooijen, 2006). Individuals with more than 50% missing genotypes and SNP loci showing evidence of strong segregation distortion (Chi-squared test; p < .00001) were excluded prior to creation of linkage maps. Maps of each linkage group (LG) were created using the Kosambi mapping function with a LOD (log of the odds) threshold of 3.0, recombination threshold of 0.40, jump threshold of 5.0, and no fixed order. A ripple was performed after each marker was added to the map. Linkage group 19, the stickleback sex chromosome (Peichel et al., 2004), was not used for the QTL analyses. All individuals and loci that were excluded from the linkage map were also excluded from QTL analyses. Final individual and marker numbers that were used in linkage map generation and QTL analyses are summarized by F1 family in Table S2.
QTL mapping was performed in R/qtl (Broman et al., 2003) on 63 traits including the 36 x- and y-coordinates from the geometric morphometrics analyses, 21 log-transformed and size-corrected linear traits from the ventral and lateral photos, log-transformed and size-corrected mean pelvic spine length (mean of left and right pelvic spine lengths) and mass, plus standard length, centroid size, right gill raker count, and total plate count (sum of left and right plate counts). QTL analyses were performed independently for each F1 family. To account for sexual dimorphism, sex was used as a covariate in all QTL analyses. Genome-wide significance thresholds were calculated for each trait using 1000 permutations, with an alpha of 0.05.
Identification of parallel QTL
For each trait, we counted the total number of unique linkage groups that contained a QTL, as well as the number of linkage groups that were found in at least two different watersheds. The latter shared linkage groups are considered “parallel” linkage groups. However, if the same linkage group contained a QTL in different F1 families of the same watershed, it was only counted once for the watershed. For each trait, we then calculated the proportion of parallel linkage groups out of the total number of unique linkage groups with a QTL. We then calculated the mean proportion of parallel linkage groups across the traits for which at least one QTL was identified in at least two different watersheds (“among-watershed parallelism”; Figure 1B). This approach will overestimate the number of parallel QTL because we were unable to directly compare linkage maps among watersheds due to the low numbers of shared markers. Therefore, we cannot exclude the possibility that the QTL map to different regions of the same linkage group. Nonetheless, given the low levels of parallelism we identified in this study, this bias is unlikely to affect our conclusions.
Using the same approach as above, we also calculated the proportion of parallel linkage groups for each pairwise comparison of watersheds (n = 6 pairwise comparisons among 4 watersheds; “between-watershed parallelism”; Figure 1B) to ask if particular pairs of watersheds had more shared QTL. Again, if the same linkage group contained a QTL in different F1 families of the same watershed, it was only counted once for the watershed. For each pairwise watershed comparison, we calculated the mean proportion of parallel linkage groups across the traits for which at least one QTL was identified in both watersheds.
We used the same approach to calculate the proportion of parallel linkage groups among F1 families of the same watershed (“within-watershed parallelism”; Figure 1B). For each among family comparison, we calculated the mean proportion of parallel linkage groups across the traits for which at least one QTL was identified in both families for Boot, Misty, and Roberts, or at least two F1 families for Pye.
To determine whether the proportion of parallel linkage groups was significantly more than expected by chance, we generated random QTL data sets, each of which was identical in structure to the real QTL data set used to calculate the observed parallelism. To generate a random QTL data set, we sampled random autosomal linkage group(s) to replace the linkage group(s) of the real data set. This was done for each trait within each F1 family. Random sampling of linkage groups was weighted by the physical lengths of linkage groups (Jones et al., 2012) to account for the fact that larger linkage groups with more genes are more likely to harbor a random QTL and that larger chromosomes have lower average crossover rates, thus biasing QTL detection to these chromosomes (Noor et al., 2001; Roesti, 2018). From each random data set, we then calculated the different types of parallelism (among-watershed, between-watershed, within-watershed) in the same way as from the empirical data. To calculate statistical significance (p values), we compared the magnitude of parallelism in 10,000 random data sets (among-watershed) or 3,000 random data sets (between-watershed, within-watershed) to the magnitude of parallelism observed in the empirical data. We also calculated the mean magnitude of parallelism across all random data sets for comparison with the observed magnitude of parallelism.
Analyses of parallel trait evolution and parallel QTL
For each of the 63 traits used for QTL mapping, we extracted the measurements of wild-caught lake and stream fish from the same watersheds (Boot, Misty, Pye, Roberts) from the Stuart et al. (2017) dataset. These included the 36 x- and y-coordinates of the 18 landmarks that match those used for QTL mapping (landmark 13 from Stuart et al. (2017) was excluded here), 21 log-transformed and size-corrected linear traits, log transformed and size-corrected mean pelvic spine length (mean of left and right pelvic spine lengths) and mass, plus log transformed standard length, centroid size, right gill raker count, left plate count, and right plate count (for this analyses, right and left plate count were considered separately). Size correction was done on the raw data from Stuart et al. (2017) following Vamosi (2002). Nomenclature of the traits in the wild-caught data was standardized to match the nomenclature used in the QTL data files. Following Stuart et al. (2017), we performed trait-by-trait linear models to test for the effects of habitat (lake vs. stream), watershed, and a habitat × watershed interaction on phenotypic divergence, with the habitat effect representing the parallel component of divergence, and the habitat × watershed effect representing the nonparallel component of divergence. Effect sizes (η2) were extracted from the linear models using the EtaSq function (BaylorEdPsych) in R.
To determine whether there is an association between the extent of phenotypic parallelism in wild fish and genetic parallelism, we calculated the Spearman’s correlation coefficient between the habitat effect size and the proportion of parallel linkage groups across all watersheds using the Hmisc package in R. Traits with either no linkage groups with a detected QTL or a single detected QTL were excluded from these analyses. For total plate count, which was mapped in the QTL analyses, the correlation was performed with the habitat effect size of left plate count only.
Results and discussion
Low genetic parallelism within, between, and among watersheds
We performed QTL mapping of 63 traits separately in each F1 family from the four lake-stream systems and found a mean of 28.6 QTL across all traits per F1 family, with a range from 6 (Roberts family 3) to 86 (Boot family 1) QTL per family (Table S2). Across all F1 families, we identified a total of 257 QTL for 56 of 63 traits (Table S3). We found evidence for both shared and unique QTL among F1 families within a watershed and across watersheds.
To quantify the extent of genetic parallelism within, between and among watersheds, we sought to identify parallel QTL. However, because we used different sets of markers for each family and watershed, it was not possible to directly compare the position of QTL across watersheds and families. Thus, we instead defined parallel QTL as those found on the same linkage group in independent watersheds or families. Although this will tend to overestimate the extent of genetic parallelism, the number of parallel linkage groups was very low for any given trait and across traits (Table 1, Table S4). Identifying parallel linkage groups allowed us to quantify the extent of genetic parallelism within, between and among watersheds for each trait, as illustrated for total plate count in Figure 1B. For example, in both Boot families (each deriving from an independent lake-stream cross from the Boot watershed), we identified a single QTL for total plate count on LG7; in this case, the proportion of parallel linkage groups within the watershed is 1. In the two Roberts families, we identified QTL on LG4 and LG17 in family 1 and a single QTL on LG4 in family 2; here, the proportion of parallel linkage groups for total plate count within the watershed is 0.5. Because Boot and Roberts have no QTL in common, the proportion of parallel linkage groups between these two watersheds is 0. Finally, of the five unique linkage groups identified across all watersheds, two (LG4 and LG7) are found in more than one watershed, leading to an among-watershed proportion of parallel linkage groups of 0.4. The same analysis was performed for each trait and then combined to examine genetic parallelism across traits among, between, and within watersheds.
Genetic parallelism between and within watersheds. For each comparison between watersheds or families within watersheds, only traits in which quantitative trait locus (QTL) were identified in both watersheds or families (or at least 2 families for the Pye comparison) were included in the analyses. Across these traits, the total number of unique linkage groups (LG) that contain a QTL in both watersheds or families, the total number of shared linkage groups (i.e., parallel) in the watersheds or families, and the mean proportion of parallel linkage groups across the traits are given. Note that the means were calculated from the proportion of parallel linkage groups calculated on a per trait basis, and therefore are not a simple proportion of the numbers given in the table. Total trait numbers are: n = 19 for the Boot–Misty comparison, n = 31 for the Boot–Pye comparison, n = 9 for the Boot–Roberts comparison, n = 19 for the Misty–Pye comparison, n = 9 for the Misty–Roberts comparison, n = 8 for the Pye–Roberts comparison, n = 18 for the Boot family 1 and 2 comparison, n = 6 for the Misty family 1 and 2 comparison, n = 15 for the Pye family 1, 2, and 4 comparison, and n = 2 for the Roberts family 3 and 5 comparison. Counts of unique and parallel linkage groups for each trait are provided in Table S4. Comparisons for which significantly more parallel linkage groups were detected by chance based on a permutation test (p < .05) are indicated with an asterisk. The full results of the permutation tests are provided in Table S5.
| Boot | Misty | Pye | Roberts | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Number unique LG | Number parallel LG | Proportion parallel LG | Number unique LG | Number parallel LG | Proportion parallel LG | Number unique LG | Number parallel LG | Proportion parallel LG | Number unique LG | Number parallel LG | Proportion parallel LG | |
| Boot | 59 | 6 | 0.180* | 79 | 10 | 0.154* | 106 | 13 | 0.154* | 30 | 5 | 0.150 |
| Misty | 13 | 2 | 0.250* | 59 | 1 | 0.026 | 27 | 3 | 0.111 | |||
| Pye | 35 | 1 | 0.033 | 22 | 1 | 0.063 | ||||||
| Roberts | 3 | 1 | 0.250 | |||||||||
| Boot | Misty | Pye | Roberts | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Number unique LG | Number parallel LG | Proportion parallel LG | Number unique LG | Number parallel LG | Proportion parallel LG | Number unique LG | Number parallel LG | Proportion parallel LG | Number unique LG | Number parallel LG | Proportion parallel LG | |
| Boot | 59 | 6 | 0.180* | 79 | 10 | 0.154* | 106 | 13 | 0.154* | 30 | 5 | 0.150 |
| Misty | 13 | 2 | 0.250* | 59 | 1 | 0.026 | 27 | 3 | 0.111 | |||
| Pye | 35 | 1 | 0.033 | 22 | 1 | 0.063 | ||||||
| Roberts | 3 | 1 | 0.250 | |||||||||
Genetic parallelism between and within watersheds. For each comparison between watersheds or families within watersheds, only traits in which quantitative trait locus (QTL) were identified in both watersheds or families (or at least 2 families for the Pye comparison) were included in the analyses. Across these traits, the total number of unique linkage groups (LG) that contain a QTL in both watersheds or families, the total number of shared linkage groups (i.e., parallel) in the watersheds or families, and the mean proportion of parallel linkage groups across the traits are given. Note that the means were calculated from the proportion of parallel linkage groups calculated on a per trait basis, and therefore are not a simple proportion of the numbers given in the table. Total trait numbers are: n = 19 for the Boot–Misty comparison, n = 31 for the Boot–Pye comparison, n = 9 for the Boot–Roberts comparison, n = 19 for the Misty–Pye comparison, n = 9 for the Misty–Roberts comparison, n = 8 for the Pye–Roberts comparison, n = 18 for the Boot family 1 and 2 comparison, n = 6 for the Misty family 1 and 2 comparison, n = 15 for the Pye family 1, 2, and 4 comparison, and n = 2 for the Roberts family 3 and 5 comparison. Counts of unique and parallel linkage groups for each trait are provided in Table S4. Comparisons for which significantly more parallel linkage groups were detected by chance based on a permutation test (p < .05) are indicated with an asterisk. The full results of the permutation tests are provided in Table S5.
| Boot | Misty | Pye | Roberts | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Number unique LG | Number parallel LG | Proportion parallel LG | Number unique LG | Number parallel LG | Proportion parallel LG | Number unique LG | Number parallel LG | Proportion parallel LG | Number unique LG | Number parallel LG | Proportion parallel LG | |
| Boot | 59 | 6 | 0.180* | 79 | 10 | 0.154* | 106 | 13 | 0.154* | 30 | 5 | 0.150 |
| Misty | 13 | 2 | 0.250* | 59 | 1 | 0.026 | 27 | 3 | 0.111 | |||
| Pye | 35 | 1 | 0.033 | 22 | 1 | 0.063 | ||||||
| Roberts | 3 | 1 | 0.250 | |||||||||
| Boot | Misty | Pye | Roberts | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Number unique LG | Number parallel LG | Proportion parallel LG | Number unique LG | Number parallel LG | Proportion parallel LG | Number unique LG | Number parallel LG | Proportion parallel LG | Number unique LG | Number parallel LG | Proportion parallel LG | |
| Boot | 59 | 6 | 0.180* | 79 | 10 | 0.154* | 106 | 13 | 0.154* | 30 | 5 | 0.150 |
| Misty | 13 | 2 | 0.250* | 59 | 1 | 0.026 | 27 | 3 | 0.111 | |||
| Pye | 35 | 1 | 0.033 | 22 | 1 | 0.063 | ||||||
| Roberts | 3 | 1 | 0.250 | |||||||||
Among-watershed parallelism
Across all watersheds and traits for which more than one QTL was identified, the mean proportion of parallel linkage groups was 0.154 (range 0 to 0.6), which was more than expected by chance (0.095, permutation P = 0.0186; Figure S3). The trait with the highest proportion of parallel linkage groups (0.6) was snout length, for which five different linkage groups (LG1, LG3, LG5, LG8, LG9) contained QTL. Three of these five linkage groups were shared across watersheds: LG1 was identified in Boot family 1 and Pye family 4, LG8 was identified in Boot family 1 and Roberts family 5, and LG9 was identified in Boot family 2 and Pye family 1 (Table S4).
Between-watershed parallelism
The proportion of parallel linkage groups varied among pairwise combinations of watersheds. The highest mean proportion of parallel linkage groups was found between Boot and the other watersheds (Boot–Misty mean = 0.154 (range 0 to 0.5); Boot–Pye mean = 0.154 (range 0 to 1); Boot–Roberts mean = 0.150 (range 0 to 0.5)), but this was more than expected by chance only in two of the three comparisons (Table 1, Table S4, Table S5). None of the other three comparisons showed a higher proportion of parallel linkage groups than expected by chance, and the lowest mean proportion of parallel linkage groups was found between Misty and Pye (mean = 0.026, range 0 to 0.5) (Table 1, Table S4, Table S5).
Within-watershed parallelism
We also found relatively low levels of parallel linkage groups among families of the same watershed (means range from 0.033 to 0.250); in some cases, this was even lower than the proportion of parallel linkage groups between watersheds (Table 1). However, the within-watershed parallelism was still more than expected by chance in two of four watersheds (Table 1, Table S5). Thus, we found that different QTL are segregating for the same trait within watersheds. To the best of our knowledge, this is one of few studies in which multiple individuals from the same population have been used to set up independent crosses in a QTL study in any system (but see Whiting et al., 2022). Two other stickleback studies have revealed that distinct QTL on LG21 have relatively large effect on branchial bone length and tooth number across multiple populations but that these QTL are not fixed within populations (Cleves et al., 2018; Erickson et al., 2018). In general, more studies that perform QTL mapping in replicate crosses from the same populations are needed to reveal the extent to which different QTL are segregating for the same trait within a given population comparison. Such genotypic redundancy has been suggested to facilitate local adaptation, particularly under high levels of gene flow (Láruson et al., 2020) and/or when traits are highly polygenic (Barghi et al., 2020; Goldstein and Holsiger, 1992; Yeaman, 2015; Yeaman et al., 2018). However, the number of redundant genotypes segregating within a population for any given phenotype (so-called segregating redundancy; Láruson et al., 2020) is mostly unknown (but see Barghi et al., 2019).
No correlation between phenotypic and genetic parallelism
Despite the relatively low proportion of linkage groups containing a QTL in more than one watershed, we investigated whether traits that showed more parallel phenotypic divergence would share more QTL. To identify the traits evolving in parallel across the four lake-stream systems used for QTL mapping, we performed linear models to quantify the consistent effects of habitat (parallel divergence) and the interaction of habitat and watershed (nonparallel divergence) on each trait measured in wild-caught lake and stream fish from these four watersheds (Stuart et al., 2017). Although 78% of the traits showed evidence for parallel divergence, with a significant effect (p < .05) of habitat in the linear model, all but three of these traits also showed a significant habitat × watershed effect (Table S4). That is, the effect on the phenotype of being from a lake or a stream depended upon the watershed of origin. Still, for 27 traits, the effect size of habitat was greater than the effect size of the interaction between habitat and watershed, suggesting parallel divergence in these traits (points above the diagonal line in Figure 2A). The two traits with the largest effect of habitat, body depth (habitat η2 = 0.497, habitat × watershed η2 = 0.024) and gill raker number (habitat η2 = 0.493, habitat × watershed η2 = 0.105), clearly showed consistent patterns of divergence between lake and stream fish across the four watersheds (Figure 2B,C). Indeed, these traits are also among the most parallel between lake-stream pairs on Vancouver Island and across the world-wide distribution of lake-stream pairs (Paccard et al., 2020; Stuart et al., 2017).
Quantification of parallel traits. (A) For each trait, the habitat × watershed (nonparallel divergence) effect size (η2) is plotted relative to the habitat (parallel divergence) effect size (η2) calculated from a linear model (Table S4). The trait means of each lake and stream population are plotted, with a line connecting the lake and stream population of each of the four watersheds, for (B) size-corrected body depth, (C) log right gill raker number, (D) size-corrected mean pelvic spine length (average of right and left pelvic spines), and (E) size-corrected snout length. This figure is a replotting of a subset of the data from Stuart et al. (2017), only including the four lake-stream pairs studied here.
Despite the phenotypic parallelism in lake-stream divergence for body depth and gill raker number, we did not find strong evidence for underlying genetic parallelism. For body depth, QTL were identified on five linkage groups, but only one of these linkage groups (LG12) contains QTL in two watersheds (Boot and Pye; Table S3, Table S4). Similarly, for gill raker number, QTL were identified on three linkage groups, but only one of these linkage groups (LG1) contains QTL in two watersheds (Boot and Roberts; Table S3, Table S4). The proportion of linkage groups that contain a QTL in more than one watershed is similar (0.333) for mean pelvic spine length, which shows the largest habitat × watershed effect size (η2 = 0.372; Figure 2D; Table S4). The trait with the highest proportion of shared linkage groups (0.60) is snout length (Table S4), which has significant (p < .05) but moderate effect sizes for both the parallel and nonparallel components of divergence (habitat η2 = 0.013, habitat × watershed η2 = 0.076; Figure 2E; Table S4). Given these few examples, it is not surprising that the proportion of linkage groups that contain a QTL for the same trait in more than one watershed is not correlated with the habitat effect size estimated from the linear models (Spearman’s correlation coefficient = 0.091; p = .561; Figure 3).
No correlation between phenotypic parallelism and genetic parallelism. For each trait, the proportion of shared linkage groups (LGs); i.e., that contain a QTL for the same trait in more than one watershed, is plotted relative to the habitat effect size (η2) for that trait, estimated from linear models. The correlation is not significant (Spearman’s correlation coefficient = 0.091; p = .561).
Caveats of our analyses
The first caveat of our analyses is that both the number of markers and the number of F2s varied among crosses (Table S2), despite our efforts to produce F2 families of similar sizes. Both the number of markers and the number of individuals impact the power to detect QTL of small effect (Beavis, 1998; Broman and Sen, 2009). However, marker numbers were relatively high (range of 230–715) and did not appear to limit power to detect QTL in our crosses, as the correlation between marker number and detected QTL was not significant (Spearman’s correlation coefficient = −0.0837, p = .831). By contrast, there was a strong correlation between the number of F2s and the number of QTL detected (Spearman’s correlation coefficient = 0.867, p = .0025). Despite this overall strong correlation, we had very similar numbers of F2s in the Boot and Misty crosses (Boot 1 = 259, Boot 2 = 274, Misty 1 = 198, Misty 2 = 214) but detected variable numbers of QTL (86, 37, 19, and 32 QTL, respectively) in these crosses (Table S2). These differences in the numbers of F2s between the crosses did not seem to impact the major conclusions of our study. For example, the highest number of QTL detected for the same trait on the same linkage group was found between crosses from the Boot and Pye watershed, although the Pye crosses had between 26% and 64% of the number of F2s in the Boot crosses (Table 1, Table S2). Furthermore, in any given comparison the overall proportion of parallel linkage groups is not correlated with the difference in the number of F2s between the crosses in that comparison (Spearman’s correlation coefficient = −0.267, p = .401; Table S6). Still, it is likely that the small numbers of F2s in some of our crosses limited our ability to detect parallel QTL (Beavis, 1998; Broman and Sen, 2009). Thus, we cannot be sure whether the lack of a QTL in a cross is because the QTL is not there or because we could not detect it. To avoid a downward bias in our estimates of genetic parallelism, we therefore included only traits for which QTL were detected in at least two families within a watershed (for within-watershed parallelism) or in at least two watersheds (for between- and among-watershed parallelism). Still, the results we present here are likely underestimates of genetic parallelism.
The second caveat of our analyses is that we needed to make separate linkage maps for each F1 family because there were very few markers that had fixed differences between the lake and stream grandparents of all the crosses for a given watershed. For example, for the Pye watershed, we only identified 62 loci that had the same homozygous SNP genotype in all three lake grandmothers and the alternative homozygous SNP genotype in all three stream grandfathers used to establish the three different Pye F1 families (Table S1), which was an insufficient number for linkage mapping and QTL analyses. Similarly low numbers of informative markers were found between the lake and stream grandparents for the crosses in the three other watersheds, which is why we required that markers only be informative within a single F1 family (Table S2). The lack of fixed SNPs between lake and stream stickleback within a watershed is consistent with the relatively low divergence between these pairs (mean genome-wide Fst: Boot = 0.065; Misty = 0.076; Pye = 0.089; Roberts = 0.072; Rennison and Peichel, 2022). Given the low number of shared fixed SNPs within a watershed and the relatively low levels of genomic parallelism in these lake-stream systems (Rennison et al., 2019), it was difficult to identify markers that could anchor the linkage maps for a direct comparison across the four watersheds. For this reason, our analysis of parallel QTL merely asks whether QTL that underlie the same trait are found on the same or different linkage group in the different watersheds. This should lead to an overestimation of the number of parallel QTL, if what we term parallel QTL are in different regions of the same linkage group. However, given the overall very low levels of parallel QTL identified in this study and the detection issues discussed in the paragraph above that would have an opposite effect on our estimates of genetic parallelism, it is unlikely that this possible bias has affected our overall conclusions.
Comparison with other stickleback systems
Despite these caveats, our results provide a striking contrast to previous results in sticklebacks, which found more repeated use of the same QTL underlying divergence in the same traits in crosses between benthic and limnetic ecotypes from two independent lakes (Conte et al., 2015) and in crosses between three independent benthic populations with the same marine population (Erickson et al., 2016). In the benthic-limnetic crosses, nearly half of the QTL identified affected the same trait in both crosses (Conte et al., 2015). Among the three benthic-marine crosses, there was a significant overlap in shared QTL, and between 40% and 47% of QTL were shared between any two crosses (Erickson et al., 2016). By contrast, only a maximum of 15% of QTL were shared between any two watersheds in our analyses, although this was greater than expected by chance in at least some of our comparisons. There are several possible explanations for the difference in the extent of genetic parallelism between benthic-limnetic and lake-stream sticklebacks, including differences in: (1) study design; (2) the source of standing genetic variation and/or evolutionary history of the populations; (3) the extent of gene flow between the ecotypes; and/or (4) the strength and/or direction of selection. We will discuss each of these possibilities in turn.
Differences in study design
First, a major difference in study design between our lake-stream study and one of the two previous studies in the benthic-limnetic system is that the Conte et al. (2015) study focused on highly parallel traits, while we studied traits chosen without prior information about the degree of parallelism. However, the Erickson et al. (2016) study did not explicitly focus on parallel traits and still observed relatively high levels of genetic parallelism. Second, we had fewer F2s per F1 family in our study relative to the Conte et al. (2015) study (n = 403 Paxton F2s; n = 323 Priest F2s), but the number of F2s (n = 186 ENOB F2s, n = 186 PAXB F2s, n = 180 PRIB F2s) in the Erickson et al. (2016) study was similar to the number of F2s in our larger F1 families (Table S2). Third, both previous studies used a single shared linkage map for QTL analyses; however, the same QTL were detected when separate linkage maps for each cross were used in the Conte et al. (2015) study. Given that similar levels of genetic parallelism were detected in the two benthic-limnetic studies despite the differences in study design between them, and the similarities in the Erickson et al. (2016) study to ours, we do not think that differences in study design can explain the relatively low levels of genetic parallelism found in our study.
Differences in the sources of standing genetic variation and/or the evolutionary history of the populations
A difference in the source of standing genetic variation in ancestral marine populations appears to contribute to low levels of parallelism between sticklebacks from the Eastern Pacific and Atlantic basins (Fang et al., 2020; Jones et al., 2012; Magalhaes et al., 2021; Roberts Kingman et al., 2021). It is unlikely that the ancestral marine source of standing genetic variation differs significantly among the benthic-limnetic and lake-stream pairs studied here, as these post-glacial freshwater populations are found within close geographical proximity on Vancouver and Texada islands within the Strait of Georgia (British Columbia, Canada). However, the evolutionary history of the benthic-limnetic systems and the lake-stream systems may differ. The benthic-limnetic species pairs found within individuals lakes likely involved a double invasion scenario, in which the lakes were colonized first by a marine ancestor that adapted to the benthic niche, followed by a second colonization approximately 2000 years later of additional marine fish, which adapted to the limnetic habitat (Taylor and McPhail, 2000). By contrast, it is not clear exactly how the colonization of lakes and parapatric streams occurred after the end of the last ice age. It is likely that they were colonized from the same ancestral population, as the lake-stream pairs studied here (with the exception of Pye) are sister taxa in a phylogenetic analysis of sticklebacks from Vancouver Island (Stuart et al., 2017). Divergence times between the lake and stream systems in this study are estimated to be 2750 (Boot), 3876 (Misty), 5729 (Pye), and 3267 (Roberts) years (Stuart et al., 2017). This is consistent with the low levels of genetic differentiation observed between lake and stream fish (mean genome-wide Fst: Boot = 0.065; Misty = 0.076; Pye = 0.089; Roberts = 0.072; Rennison and Peichel, 2022). Genetic differentiation is much higher between benthic and limnetic fish (mean genome-wide Fst: Paxton = 0.206; Priest = 0.197; Rennison and Peichel, 2022), which may suggest an older divergence time (or stronger divergent selection; see below) between benthics and limnetics. Older divergence times would be consistent with a scenario in which sticklebacks had previously adapted to benthic or limnetic freshwater habitats and old adaptive alleles were maintained as standing variation in the marine populations that founded the extant benthic-limnetic species pair lakes, as has been observed for alleles that contribute to marine-freshwater divergence (Nelson and Cresko, 2018; Roberts Kingman et al., 2021). If benthic-limnetic divergence is due to selection on such ancient standing genetic variation, we would expect the higher level of genetic parallelism observed in this system (Conte et al., 2015; Erickson et al., 2016).
Differences in the extent of gene flow between the populations
Higher levels of gene flow might select for a clustered genetic architecture in which adaptive loci are linked in regions of low recombination, thereby constraining the location of QTL and leading to higher levels of parallelism (Feder et al., 2012; Kirkpatrick and Barton, 2006; Yeaman, 2013; Yeaman and Whitlock, 2011). In addition, there may be a bias for detecting QTL in regions of low recombination (Noor et al., 2001; Roesti, 2018). Clustering of QTL in regions of low recombination, particularly on chromosomes 4, 7, and 21, has been observed across all stickleback QTL data, as well as within individual studies using crosses derived from benthic populations (Arnegard et al., 2014; Erickson et al., 2016; Miller et al., 2014; Peichel and Marques, 2017). We observe less clustering in our data, although linkage group 4 has more QTL than expected given either the number of genes or the length of the chromosome (Table S7). Some evidence of a clustered genetic architecture is also suggested by the observation that a greater number of unique traits map to linkage groups with a higher proportion of parallel QTL (Spearman’s correlation coefficient = 0.583, p = .007; Table S8). These findings further imply that parallel QTL are more pleiotropic (either due to linkage of multiple mutations or due to a single pleiotropic mutation), consistent with the conclusions of a population genomic study in the same lake-stream systems, which found that genomic regions differentiated in multiple lake-stream pairs contained more QTL than those differentiated in a single lake-stream pair (Rennison and Peichel, 2022). However, this same study also found more repeatable and clustered patterns of genomic divergence in the benthic-limnetic systems than in the lake-stream systems. Overall, these differences in the extent of genomic clustering and genetic parallelism are consistent with the extent of gene flow in the two systems: gene flow is on the order of 10−4 to 10−5 individuals per generation in the lake-stream systems (based on whole-genome RAD-seq data; Stuart et al., 2017), and an order of magnitude higher in the benthic-limnetic systems (based on a handful of microsatellite markers; Gow et al., 2006). Although these gene flow estimates rely on different data types and analytical methods, they suggest that differences in gene flow might contribute to differences in the extent of genetic (and genomic) parallelism between the systems.
Differences in the strength and/or direction of selection
Theoretical work suggests that the probability of parallel genetic evolution is higher under natural selection than under neutrality (Orr, 2005). A more recent modelling study demonstrated that genetic parallelism is highest when selection is completely parallel but decreases rapidly as selection becomes less parallel, particularly when adaptation occurs from standing genetic variation (Thompson et al., 2019). Even when selection is completely parallel, there is a reduction in genetic parallelism when populations start at a greater distance from the phenotypic optimum than when they start closer to the optimum (Thompson et al., 2019). Interestingly, both the direction and magnitude of phenotypic parallelism is much higher in benthic-limnetic than in lake-stream pairs of sticklebacks in British Columbia (Oke et al., 2017), suggesting that divergent selection might be more parallel in the benthic-limnetic than in the lake-stream systems. The higher Fst, despite higher levels of gene flow, observed in the benthic-limnetic systems also suggests that divergent selection is stronger in the benthic-limnetic than in the lake-stream systems. Although this hypothesis requires future empirical tests, we currently favor the possibility that differences in the direction and magnitude of selection between the benthic-limnetic and lake-stream systems, perhaps coupled with differences in evolutionary history and gene flow, might explain the reduction of genetic parallelism in the lake-stream stickleback systems.
Comparison with other systems
Although genetic parallelism seems to be prevalent when repeated phenotypic evolution is due to loci of large effect (Conte et al., 2012; Martin and Orgogozo, 2013; Stern, 2013), it has been predicted that lower levels of genetic parallelism should be found for the repeated evolution of phenotypes with a polygenic architecture (Orr, 2005; MacPherson and Nuismer, 2017). However, surprisingly few comparative QTL mapping studies of repeatedly evolved traits with an underlying polygenic architecture have been done in any system. In addition to the studies in threespine stickleback (Conte et al., 2015; Erickson et al., 2016; this study), comparative QTL mapping studies have examined the genetic architecture of repeated pelvic reduction in three populations of ninespine stickleback (Kemppainen et al., 2021), repeated evolution of leaf shape in three closely-related Mimulus species (Ferris et al., 2015), repeated divergence of male mating song in three species pairs of Laupala crickets (Blankers et al., 2019), and repeated evolution of mimetic wing patterns in two species of Heliconius butterflies (Bainbridge et al., 2020). The findings across these disparate studies are mixed. As in the two benthic-limnetic stickleback studies, more genetic parallelism than expected by chance was found in the Mimulus and Laupala studies (Blankers et al., 2019; Ferris et al., 2015). As in the current lake-stream stickleback study, relatively low levels of genetic parallelism are seen in the ninespine stickleback and Heliconius studies (Bainbridge et al., 2020; Kemppanien et al., 2021). Interestingly, the relatively low levels of parallelism observed for small effect loci in the Heliconius study contrasts to the reuse of a major effect locus (optix) for color patterns in these same species (Reed et al., 2011), consistent with the prediction that loci of large effect should be more parallel than those of small effect (MacPherson and Nuismer, 2017; Orr et al., 2005). Although, we could not test this prediction in our study because we only identified loci of relatively small effect, two other studies have found no correlation between QTL effect size and genetic parallelism (Blankers et al., 2019; Conte et al., 2015). However, more comparative QTL mapping studies of repeatedly evolved phenotypes with both simple and polygenic architectures are required to test whether genetic parallelism is higher for loci of large effect.
Conclusions
The threespine stickleback is well-known for remarkable cases of parallel evolution at the genetic level, facilitated by repeated selection on standing variation (e.g., Colosimo et al., 2005; Howes et al., 2017; Indjeian et al., 2016; Kitano et al., 2010; Miller et al., 2007) or by recurrent mutation at the same locus (Chan et al., 2010; Shapiro et al., 2004; Xie et al., 2020). However, recent studies of specific traits have also revealed that different loci can underlie the repeated evolution of similar phenotypes in threespine stickleback (Ellis et al., 2015; Glazer et al., 2015). Our more comprehensive analyses of many traits in QTL crosses from four different lake-stream pairs emphasize that repeated evolution at the genetic level is not the only mechanism that underlies repeated phenotypic divergence in sticklebacks. Furthermore, our results show that there is no correlation between genetic and phenotypic parallelism, suggesting that genetic constraints are unlikely to play a major role in explaining phenotypic parallelism and that there is a high level of genotypic redundancy even within stickleback populations. Although the reasons for differences among stickleback systems or traits are unknown, our study suggests that future research should move from simply documenting patterns of parallel evolution at the genetic and genomic level to identifying the drivers of these patterns.
Data availability
All datasets and scripts used for analyses are deposited on Dryad (doi:10.5061/dryad.63xsj3v5r). Raw double-digest RAD sequencing reads are available from the NIH Sequence Read Archive (BioProject PRJNA817082).
Author contributions
A.P.H., D.I.B., and C.L.P. conceived and designed the study. C.L.P. and A.P.H. collected the fish for the crosses. C.L.P. made the crosses. H.A.P. phenotyped the F2s. H.A.P. and Y.E.S. made the ddRADseq libraries with input from D.I.B. H.A.P., D.J.R., M.R., and C.L.P. analyzed the data. C.L.P. wrote the paper with input from all authors.
Funding statement
The project was supported by US National Science Foundation grants DEB-1144773 (D.I.B., A.P.H.) and DEB-1144556 (C.L.P.).
Conflict of interest: The authors declare no conflict of interest.
Acknowledgments
We thank Rowan Barrett, Matt Dubin, Susannah Halbrook, Dieta Hanson, and Carole Tanner for help in the field, Shaugnessy McCann for fish care, Andy Marty of the Fred Hutchinson Cancer Research Center Genomics Shared Resource for assistance with Illumina sequencing and data transfer, and Zuyao Liu for help with upload to the NIH Sequence Read Archive.
