Abstract
Heritable phenotypic differences between populations, caused by the selective effects of distinct environmental conditions, are of commonplace occurrence in nature. However, the actual genomic targets of this kind of selection are still poorly understood. We conducted a quantitative trait locus (QTL) mapping study to identify genomic regions responsible for morphometric differentiation between genetically and phenotypically divergent marine and freshwater nine-spined stickleback (Pungitius pungitius) populations. Using a dense panel of SNP-markers obtained by restriction site associated DNA sequencing of an F2 recombinant cross, we found 22 QTL that explained 3.5–12.9% of phenotypic variance in the traits under investigation. We detected one fairly large-effect (PVE = 9.6%) QTL for caudal peduncle length–a trait with a well-established adaptive function showing clear differentiation among marine and freshwater populations. We also identified two large-effect QTL for lateral plate numbers, which are different from the lateral plate QTL reported in earlier studies of this and related species. Hence, apart from identifying several large-effect QTL in shape traits showing adaptive differentiation in response to different environmental conditions, the results suggest intra- and interspecific heterogeneity in the genomic basis of lateral plate number variation.
Similar content being viewed by others
Introduction
Adaptation to different environmental conditions is often, but not always1, accompanied by genetically based morphological divergence in size and shape2,3,4. While common garden experiments5,6,7 can verify the heritable nature of such divergence, uncovering the genetic basis of these complex phenotypic traits can be far more challenging8,9. For instance, adaptive genetic divergence in body shape among fish populations residing in different environments has been repeatedly demonstrated10,11,12, but the genetic underpinnings of this divergence are still fairly poorly understood13,14. This is not surprising, because body shape is a complex trait, likely to be highly polygenic: large sample sizes, both in terms of number of individuals and markers, are needed to identify the causal loci influencing variation in such traits15,16. The quest for understanding the evolution of body shape is further complicated by the fact that different aspects of body shape variability can be under conflicting selection pressures, and genetic correlations caused by pleiotropy and linkage disequilibrium can constrain or facilitate allele frequency changes in a given locus depending on the prevailing selection presures17.
The stickleback fishes (Gasterosteidae) provide excellent model systems for studies of the genetic architecture of body shape divergence. The three-spined stickleback (Gasterosteus aculeatus) has in fact been proposed as a model to study the evolution of body shape in fish17; recently the nine-spined stickleback (Pungitius pungitius)–which diverged from the three-spined stickleback around 13 million years ago18–has also been emerging as a model for evolutionary investigations19. These two species are ecologically, and to a certain degree also phenotypically, very similar20,21. Early mapping studies in three-spined sticklebacks have focused on simple morphological traits, such as pelvic reduction22,23 and armor loss24,25,26, and have been followed by studies focusing on the genetic architecture of complex traits including body shape variability13,27,28,29,30. Similarly to three-spined sticklebacks, freshwater populations of nine-spined sticklebacks have repeatedly and independently evolved deeper bodies, reduced armor, shorter caudal peduncles, smaller brains, and different behavioral characteristics as compared to marine populations31,32,33. Quantitative genetic studies conducted in ‘common garden’ conditions suggest an additive genetic basis for these morphological divergences7,22,34,35,36,37. Thus, nine-spined sticklebacks not only provide another promising model to examine the genetics of body shape divergence as an adaptation to life in freshwater environments, but also offer a chance to explore whether the genetic basis (i.e. common genes and/or genetic pathways) of body shape divergence is similar to that in the three-spined stickleback. Addressing this question can provide important insights into the potential role of genetic and developmental constraints in the evolution of complex phenotypes.
Identifying the genomic regions that control phenotypic variation is the first step towards understanding the genetic underpinnings of adaptive divergence among populations38,39. Quantitative trait locus (QTL) mapping is a classical method used for this aim. In the past, QTL-mapping in non-model organisms has relied on low-density genetic maps, typically comprised of only a few hundred molecular markers. Advances in genomic technologies have made it feasible to explore the genetic architecture of phenotypic traits at a genome-wide scale in both model and non-model species40,41,42,43. There are now several approaches (e.g., multiplexed shotgun genotyping44, reduced-representation sequencing45,46, and restriction-site associated DNA sequencing [RAD-seq]47,48) that allow the discovery and genotyping of thousands of markers across any genome of interest, even in non-model organisms with limited or no genomic resources48,49. RAD-seq has been utilized to construct high-density linkage maps and to detect QTL in an increasing number of studies14,50,51,52. Although the power and precision of QTL-mapping critically depends on the experimental design and number of mapped progeny53, Stange et al.54 demonstrated that high-density maps can increase the precision of QTL localization and effect sizes, especially for small and medium sized QTL, as well as the power to resolve closely linked QTL.
The main goal of this study was to investigate the genetic architecture of morphometric divergence between marine and freshwater nine-spined sticklebacks, with the aid of QTL-mapping using thousands of SNPs obtained with a RAD-seq approach. To this end, we used 283 F2-generation full-sib offspring derived from a F1-generation inter-cross between phenotypically and genetically divergent marine and freshwater populations. We mapped genomic positions in a total of 49 traits, including 38 principal component (PC) scores for body shape, 10 anatomical morphometric traits (Fig. 1), as well as lateral plate number. The detected QTL were compared to those observed in earlier studies of sticklebacks to see whether the same or different QTL for homologous traits were discovered in different studies. We also conducted functional annotation of the QTL regions in order to identify candidate genes controlling variation in studied traits.
Landmark positions: 1, Anterior extent of maxilla; 2, Posterior extent of supraoccipital; 3, Anterior insertion of first dorsal spine; 4, Anterior insertion of dorsal fin; 5, Posterior insertion of dorsal fin; 6, Origin of caudal fin membrane on dorsal midline; 7, Posterior extent of caudal peduncle; 8, Origin of caudal fin membrane on ventral midline; 9, Posterior insertion of anal fin; 10, Anterior insertion of anal fin; 11, Insertion point of pelvic spine into the pelvic girdle; 12, Posterior extent of ectocorocoid; 13, Anterior extent of ectocorocoid; 14, Posterior-dorsal extent of operculum; 15, Posterior-ventral extent of preopercular; 16, Dorsal extent of preopercular; 17, Posterior extent of orbit; 18, Ventral extent of orbit; 19, Anterior extent of orbit; 20, Anterior-ventral extent of preopercular; 21, Posterior extent of maxilla. Definitions of metric traits: a, head length; b, upper jaw length; c, lower jaw length; d, orbit diameter; e, dorsal fin base length; f, anal fin base length; g, caudal peduncle length; h, caudal peduncle width; i, body depth; j, snout length; k, standard body length. Measurement data and photos were collected by J. Y.
Results
Linkage map
The linkage map used in this study included 14,998 unique SNP markers distributed across 21 LGs matching the expected number of chromosomes (2n = 42;55) in the nine-spined stickleback, and was adopted from Rastas et al.56. The sex-averaged map spanned 2,529 cM, with 5.99 markers/cM56. In this map, 5,241 SNPs distributed on 4,791 reference sequences could be uniquely mapped to the three-spined stickleback genome, and were defined as informative SNPs56. To overcome computational limitations in QTL mapping due to high-density of makers, a coarse-mapping was first conducted with a simplified linkage map with 466 informative SNPs (Supplementary Table 1), which was followed by fine-mapping with additional SNPs around QTL regions identified by the coarse-mapping (see Methods for details).
Morphological variation
The original measurements of morphological traits in the F2 progeny are given in Supplementary Table 2. Principal Component Analysis (PCA), based on landmark positions, was used to identify the independent axes of body shape variation. This analysis identified 38 PCs, of which the first three each accounted for >10% of the total body shape variation (Supplementary Fig. 1). PC1, accounting for 35.2% of the total variation, captured primarily variation in body depth and caudal peduncle length. Along this axis, the F2 progeny varied from individuals having shallow bodies and long caudal peduncles to individuals having deep bodies and short caudal peduncles (Fig. 2). PC2 and PC3, accounting for 15.7% and 10.4% of the total variation, respectively, captured variation not only in the body depth and caudal peduncle length, but also in shape variation corresponding to bending of the body downwards (PC2) and upwards (PC3; Supplementary Fig. 1). Sexual dimorphism along PC1 and PC3 was evident (t-tests, t281 ≥ 7.36, P < 0.001; Fig. 2), and thus, sex was used as a covariate to control the potential gender dependent variation in the subsequent QTL-mapping analyses. No sexual dimorphism was detected in PC2 (t-test, t281 = −0.87, P > 0.05).
Scatterplot of the first two principal component axe based on analysis of all landmarks. Black dots depict males, and red dots depict females. Wireframe graphs illustrate the body shape variation along the first principal component axis; black dots in the wireframes indicate the 21 landmarks used in shape analyses.
The ten continuous traits that were used in the QTL-mapping showed obvious divergence between wild collected marine (HEL) and pond (RYT) fish (Supplementary Fig. 2). For example, marine sticklebacks had longer caudal peduncles, narrower bodies, shorter lower jaws and snouts than pond individuals (Supplementary Fig. 2). Likewise, marine individuals had on average more lateral plates than pond individuals (Supplementary Table 3). Sexual dimorphism was evident in most of the traits (results not shown;31), and thus sex was included as covariate to all QTL-analyses.
QTL-mapping
With coarse-mapping we detected a total of 22 QTL (ten for PC scores of body shape, four for anatomical measures, and eight for lateral plate number) on 11 different LGs, that were significant at the genome-wide level (Supplementary Table 4). The significant QTL for anatomical measures were associated with lower jaw length, caudal peduncle length, body depth, and snout length, whereas no significant QTL were found for the six other measures (Supplementary Table 4). After adding more markers around QTL regions detected with the simplified map, the linkage map for fine-mapping was significantly improved in terms of marker density–approximately 2.1 markers/cM. With fine-mapping, more accurate (as judged from narrower CIs) marker positions for each of the 22 QTL regions were obtained (Table 1; Figs 3, 4, 5). Eighteen of the 22 significant QTL identified in the coarse-mapping became replaced by new and more accurate QTL makers in the fine-mapping results (Table 1; Supplementary Table 4).
Significant QTL are marked by different colors. The QTL bars represent 1.5 unit confidence intervals. The graphs on the right side (Y-axis) of each linkage group show LOD score distribution, with dotted threshold line.
Significant QTL are marked by red color. The QTL bars represent 1.5 unit confidence intervals. The graphs on the right side (Y-axis) of each linkage group show LOD score distribution, with dotted threshold line.
Significant QTL for left, right and total lateral plate numbers are marked by different colors. The QTL bars represent 1.5 unit confidence intervals. The graphs on the right side (Y-axis) of each linkage group show LOD score distribution, with dotted threshold line.
QTL for body shape variation
While ten QTL markers on seven LGs showed significant association with nine PC scores of body shape variation in the coarse-mapping (Supplementary Table 4), seven of these markers were refined in the fine-mapping analyses (Table 1; Fig. 3). The percentage of variance explained (PVE) by the individual QTL varied from 3.50% to 12.90% (Table 1). A large effect QTL (PVE > 10%) was detected on LG7 (6.98 cM) for two PC scores (PC6 and PC11). In addition to the large effect QTL, another QTL was found on LG7 for PC3 (15.48 cM), which was also affected by a QTL on LG8 (76.07 cM). Two QTL on LG17 were associated with PC16 and PC20 (31.00 and 43.35 cM, respectively). QTL were also found on LG4 (78.06 cM), LG5 (27.59 cM), LG14 (87.17 cM), and LG15 (22.66 cM) for PC14, PC13, PC33, and PC1, respectively (Table 1).
QTL for anatomical measures
All of the four QTL for anatomical measures detected in coarse-mapping were also retrieved by the fine-mapping, which yielded significant (at genome-wide level) QTL for body depth, caudal peduncle length, lower jaw length, and snout length (Table 1; Fig. 4). Except in the case of lower jaw length, fine-mapping refined the QTL-positions obtained from the coarse-mapping. A QTL for body depth was on LG4 (56.87 cM) with the PVE of 8.60%, a QTL for lower jaw length on LG19 (105.58 cM) with PVE of 6.70%, a QTL for snout length on LG20 (46.8 cM) with PVE of 7.40%, a QTL for caudal peduncle length on LG15 (12.11 cM) with the PVE of 9.60% (Table 1). The F2 progeny with different QTL genotypes in these loci differed significantly in their mean phenotypic values (Fig. 6): individuals with AA genotype (pond allele) on marker 27323 had significantly longer lower jaws (ANOVA-LSD, P < 0.01) than CC-homozygotes (marine allele); on marker 13320, individuals with CC genotype (pond allele) had significantly shorter caudal peduncles (ANOVA-LSD, P < 0.01) than GG-homozygotes (marine allele; Fig. 6). Likewise, individuals with AA genotype (pond allele) on marker 11319 had significantly deeper bodies (ANOVA-LSD, P < 0.01) than GG-homozygotes (marine allele); on marker 21583, individuals with CC genotype (pond allele) had significantly longer snouts (ANOVA-LSD, P < 0.01) than TT-homozygotes (marine allele; Fig. 6; Supplementary Table 5).
In each figure, genotype in left refers to pond genotype and that on right to marine genotype. Box indicates the lower and upper quartile values and whiskers represent the extreme values. Outliers are displayed by circles.
QTL for lateral plate number
At the genome-wide level, coarse mapping detected four, two, and two significant QTL for left, right, and total lateral plate number, respectively (Supplementary Table 4). In fine mapping, each of these QTL markers were replaced by new markers (Table 1; Fig. 5). QTL on LG20 (53.97 cM) and LG21 (84.89 cM) were significantly associated with variation in the left side, right side, and total plate number counts with PVEs ≥8.00% (Table 1). QTL significant at the genome-wide level were also found for the left (but not right) side plate numbers on LG8 (65.63 cM) and LG12 (76.76 cM, Table 1; Fig. 5). However, when mapping was performed for each chromosome independently, a suggestive QTL for right side and total plate number was detected on LG8 (65.63 cM), together with a suggestive QTL for total plate number on LG12 (76.76 cM). A small but evident LOD-score peak for right side plate number was observed on LG12 (76.76 cM) though it did not pass the significance threshold (Supplementary Fig. 3). Genotypes of the marker 12832 for the QTL on LG8 were missing in both of the grandparents, and thus origin of the alleles could not be identified. The F2 progeny with different QTL genotypes in the remaining three loci differed in their mean plate number values (Fig. 7): individuals with CC genotype (pond allele) on marker 22134 had fewer plates (ANOVA-LSD, left plate number: P < 0.05; right plate number: P > 0.05; total plate number: P < 0.05) than TT-homozygotes (marine allele); on marker 11482, individuals with TT genotype (pond allele) had significantly fewer plates (ANOVA-LSD, left plate number, right plate number, and total plate number: P < 0.01) than CC-homozygotes (marine allele); on marker 18769, individuals with CC genotype (pond allele) had significantly fewer plates (ANOVA-LSD, left plate number, right plate number, and total plate number: P < 0.01) than AA-homozygotes (marine allele; Fig. 7). Although several QTL were detected for lateral plate number variation, none of these were found on LG4, which contains the Eda-gene known to be the major gene controlling lateral plate variation in the three-spined stickleback24,57.
In each figure, genotype in left refers to pond genotype and that on right to marine genotype. Box indicates the lower and upper quartile values and whiskers represent the extreme values. Outliers are displayed by circles.
Candidate genes for morphological divergence
All of the QTL regions identified in the fine-mapping could be annotated with the aid of the three-spined stickleback genome (Table 1). These regions contained several genes with functions related to, for example, RNA or protein binding, and cellular protein catabolic processes (Table 1). Most of the detected QTL markers were located outside of these genes; only one QTL for PC6 and PC11 was located within the Ccdc90b-gene (Table 1). In total, 92 genes were identified within the (1.5) confidence interval for a given QTL for shape variation, and 49 genes were identified for lateral plate number variation (Table 1). Pathway identification with the KEGG database showed that these 49 genes participate in several pathways, but none of them is involved in the Eda-pathway (cytokine-cytokine receptor interaction pathway) controlling for lateral plate number variation in the three-spined stickleback24,57.
Discussion
The most important findings of this study include identification of several fairly large-effect QTL for phenotypic traits of ecological, evolutionary and systematic importance. Hence, the results can advance our understanding of the genetics and evolution of traits of major adaptive and systematic significance. Furthermore, by identifying three novel QTL (located on LG8, LG20 and LG21 respectively) associated with variation in lateral plate numbers which differ from those detected in earlier study of this58 and related species23,24, the results suggest a heterogeneous genomic basis for a trait of major evolutionary and systematic significance. Lateral plate number is an important diagnostic trait in taxonomy and systematics in the genus Pungitius59,60, and hence, our findings are relevant for determining the utility of this trait for systematic inference. Apart from these findings and considerations, the results serve to illustrate the possibilities and challenges associated with QTL-mapping with large number of markers as generated by RAD-seq. In the following paragraphs, we will discuss each of the above mentioned points in light of our findings and related issues.
An earlier study of morphometric divergence between pond and marine nine-spined sticklebacks–including the two populations used in this study–has suggested this divergence has a genetic basis31. Among other things, marine nine-spined sticklebacks have narrower and longer caudal peduncles than the pond fish. The caudal peduncle is associated with maneuvering and locomotion performance in fish61,62,63. The elongate peduncle increases the amplitudes needed to drive the caudal fin and allows for the control over the angle of attack of the caudal fin64,65, and is likely adaptive for fish moving in open water environments and under high predation risk. There is also considerable activity in the caudal peduncle when a fish changes its direction of movement65. Here, we identified a significant QTL on LG15 for the variation of peduncle length between marine and pond populations, indicating that the allele from marine population contributes to elongation of the caudal peduncle. The results suggest that the detected QTL may provide a starting point to decipher the genetic underpinnings and molecular mechanisms of adaptive divergence in caudal peduncle length in sticklebacks. Given that the strength of predator-mediated selection on caudal peduncle length can be quantified easily in mesocosm settings in sticklebacks66, the discovery of fairly large effect QTL in caudal peduncle length might also provide an opportunity to study dynamics of genetic variation in QTL under directional selection.
Variation in the number of lateral armor plates in stickleback fishes has received considerable attention for at least two reasons. First, being a conspicuous, variable and easily studied trait, variation in plate numbers has been used as a diagnostic trait in stickleback systematics59,60. Second, at least in the three-spined stickleback, the adaptive value of variation in lateral plate numbers is fairly well understood67, and several studies have demonstrated the adaptive nature of temporal68,69 and spatial24,70,71,72 variation in plate numbers. While the genetic basis of lateral plate number variation in the three-spined stickleback is controlled by a major QTL in the locus close to the Eda-gene, together with several minor QTL24,57, different large effect QTL have been identified to control variation in lateral plate numbers in North American nine-spined sticklebacks58. In this study, we detected two large QTL (located on LG20 and LG21) for lateral plate numbers, which were different from those discovered in North American nine-spined sticklebacks (Table 2). None of the genes located near these two QTL are involved in the Eda-pathway (cytokine-cytokine receptor interaction pathway) according to current gene annotation information, suggesting that the genetic mechanisms controlling for lateral plate variation may be even more variable than previously anticipated. To this end, our findings give support to the view that similar morphological changes might be commonly achievable through different QTL and/or genetic pathways73,74, albeit more interpopulation crosses and families would be needed to verify such a conclusion.
If the genetic basis of variation in armor traits in sticklebacks frequently differs from one population and species to another, and is subject to recurrent losses and gains over short evolutionary time scales69, lateral plate phenotypes may carry little information about systematic relationships among different taxa. Hence, the application of this trait in Pungitius systematics (see: ref. 60 for a review) may not be warranted. We also note the QTL on LG8 that influences variation in lateral plate numbers on the left side of the body was near to the QTL region influencing body shape (Figs 3 and 5). Earlier quantitative genetic30 and QTL-mapping studies13,27 have observed genetic links between shape and armor traits. Such results might help to explain why body shape differentiation in sticklebacks is often accompanied by plate number differentiation during evolutionary adaptation to freshwater environments.
In this study, we identified ten significant QTL contributing to divergence in body shape, and an additional 12 QTL contributing to variation in anatomical morphological traits and lateral plate number. All of the detected QTL had fairly large PVE values (average PVE = 8.48%) and some can be considered as large effect QTL (PVE > 10%) according to conventional standards75,76,77. A notable feature of our results is that for most traits–with the exception of PC3 and lateral plate numbers–only one single QTL was detected for each trait. While such results could be interpreted to suggest that single genes with large effects, rather than many genes with small effects, contribute to the observed phenotypic variability, such a conclusion may not be warranted from our data. Namely, the possibility that many genes with small effects contribute to the shape variation cannot be dismissed, as QTL studies are biased towards detecting QTL with large effects78,79. For instance, although we used a large number of markers, the modest size of our experiment in terms of number of F2-progeny (from a single family) may not have allowed the detection of many small effect QTL80,81. Furthermore, our decision to use stringent genome-wide significance as a criterion for calling QTL lead to the exclusion of many (n = 95) QTL which reached significance only at a chromosome-wide level. We believe that their exclusion from further considerations was justified given the statistical, and thereby also biological, uncertainty associated with them. It should also be pointed out that variation in shape is a cumulative effect of variation in multiple principal components, and hence of multiple QTL, even if variation along each individual principal component axis would be coded by a single or few QTL. Considering all these points, our results are not at odds with the view that complex morphological traits, such as shape, are likely to often have a polygenic basis82,83,84.
Previous studies have shown that different aspects of shape and morphology, such as lateral plate numbers, have evolved in similar directions in different freshwater populations of sticklebacks11,13,31,58,85. Such parallel evolution of trait complexes would not be likely if there were strong antagonistic genetic correlations among traits selected to change in a parallel fashion. However, quantitative genetic studies of sticklebacks suggest positive genetic correlations among, for instance, lateral plate numbers and several shape traits30. The ultimate source of these genetic correlations is pleiotropy and physical linkage among loci influencing variation in different traits. In this study, we found that one QTL region on LG7 (6.79–6.98 cM) affected two (by definition independent) principal component scores (PC6 and PC11; Table 1). This observation suggests that the same genes or genetic regions can control different components of shape variation, a characteristic that might facilitate rapid population divergence in shape. Likewise, a short genomic region on LG20 (46.61–53.97 cM) was associated with variance in both lateral plate numbers and snout length, indicating that the same genetic factor(s) may govern (part of) the variability in these two traits. However, whether a single pleiotropic gene or multiple linked genes control variation in both traits cannot be assessed from our data. The same applies to the interpretation of QTL for each individual PC-axis: since the shape variation captured by each PC-axis captures variance in multiple landmark coordinate positions, a QTL for a given PC-axis can be inferred to have pleiotropic effects on multiple landmark positions.
For all of the 22 QTL we detected, the precision of the QTL locations were very accurate, as judged from the narrow confidence intervals around the QTL positions. This high precision is also apparent if we compare the average width of the confidence region in this study with those of the earlier QTL studies of sticklebacks (Supplementary Fig. 4). The high precision of the QTL regions in this study is likely due to the higher density of markers than any of the earlier studies, as well as the fine-mapping approach, which narrowed the confidence intervals for the QTL positions (compare CIs in Table 1 and Supplementary Table 4). Interestingly, the high precision of QTL locations allowed us to discover that most (21/22) of the QTL were located in non-genic regions. We can identify three possible explanations for this. First, it could be that the RAD-seq method is biased towards finding polymorphic SNPs adjacent to genes rather than within. For instance, if there are more restriction sites outside than within genic regions, this would lead to a bias in detecting more non-genic than genic QTL. However, comparing the distribution of restrictions sites in the three-spined stickleback genome does not suggest such a bias: the number of PstI restriction sites on a given chromosome was significantly correlated with chromosome length (rs > 0.98, P < 1.71 × 10−16), and about 63% of the restriction sites were located in the genic regions. Second, genic regions of a genome are known to evolve under more stringent constraints than non-genic ones86,87,88,89,90, and therefore, the likelihood of detecting polymorphic SNPs outside of genes may be increased. A third and mutually nonexclusive possibility is that the variation associated with the detected QTL is not controlled by a sequence polymorphism within the genes, but in the regulatory regions outside of the genes. Our data do not allow us to disentangle these alternatives. However, given the increasing evidence for the importance of non-genic regulatory elements in controlling phenotypic variation91,92,93, it is possible that the majority of the detected QTL represent regulatory polymorphisms.
Most of the genes identified in each QTL region were classified into broad GO categories and pathways. Unfortunately, the exact function of most of these genes is poorly known. For instance, the QTL for PC6 and PC7 were mapped within the Ccdc90b gene, whose function is still uncharacterized. However, Meis homeobox 2a (Meis2a) and limb bud and heart homolog (Lbh) genes in the QTL region for PC1 are known to be primarily involved in the formation of the viscerocranium and craniofacial morphogenesis in zebrafish94 and cichlid fish95, respectively. Similarly, Fgf receptor-like 1a (Fgfrl1a) gene in the QTL region for PC14 was found to be necessary for cartilage formation in zebrafish96. The gene Fxr1 in the QTL region for PC3 is an RNA binding protein that plays a critical role in eye development and cranial cartilage derived from cranial neural crest cells97. Hence, Meis2a, Lbh, Fgfrl1a, and Fxr1 might be involved in adaptive evolution of stickleback by contributing to the regulation of cranial shape. In addition, since bone morphogenetic protein 4 (Bmp4) in the QTL region for PC1 has been reported to play an important role in the formation of the dorsal-ventral pattern in zebrafish98,99, this gene might be responsible for shape variation in the nine-spined stickleback. Another gene, identified in the QTL region for lateral plate number variation, codes for tRNA methyltransferase 1 like (Trmt1l) protein and has a wide set of functions including metal ion and RNA binding. The Trmt1l gene has been reported as a significant regulator in motor coordination and exploratory behavior in murine studies100. However, although associations between lateral plate number and behavioural variation in sticklebacks are known, it is not immediately obvious how and why variation in Trmt1l gene is associated with plate number variability in our cross. In addition, we found that the QTL associated with caudal peduncle length was located close to thioredoxin-related transmembrane protein 1 (Tmx1) gene. Unfortunately, little is known about the precise function of these two genes: Tmx1 gene is a member disulphide isomerase gene family regulating highly conserved enzyme-mediated disulphide bond formation affecting over one-third of all eukaryotic proteins101.
Most of the earlier QTL-mapping studies of sticklebacks have utilized low-density microsatellite marker-based linkage maps13,17,23,27,28,58,71,102,103,104,105,106,107. Here, we used a high-density SNP-based linkage map generated by RAD-seq technology, which allows genotyping a very large number of SNP markers for many individuals in a single step40,47,48. The RAD-seq approach has been utilized to construct linkage maps for QTL-mapping purposes in several other species38,50,51,52,108,109, including various teleost fishes14,110,111,112,113,114,115,116. However, most of these earlier studies have used a modest number of markers (range = 436–8,790; median = 2,011) compared to the linkage map used in our study. Using a high number of markers poses challenges for both linkage map construction and QTL mapping. As for the linkage map construction, many of the available software are not fully automated and require considerable user involvement in map construction (but see: e.g. ref. 56). As for QTL mapping, we used a two-stage approach (cf. coarse + fine mapping) to overcome the computational challenges associated with large marker numbers. Although this approach is not expected to improve neither the QTL detection power nor the proportion of explained phenotypic variance54, it allowed us to obtain narrow confidence intervals and hence increase the precision of QTL locations. This strategy in combination with the multiple QTL mapping approach allowed us to further filter out likely false QTL peaks detected by simple (univariate) interval mapping. In fact, the relatively low number of significant QTL reported in this study likely owes to the facts that (1) we reported only the QTL significant at the genome-wide level, and that (2) some of the QTL detected with interval mapping disappeared in multiple QTL mapping. However, given this stringent approach for QTL detection, we feel fairly confident that the QTL we reported are not only real, but also very precise.
We used Procrustes shape coordinates to map shape variability in our data. While these variables derived with geometric morphometric methods have many virtues117, they may be sensitive to distortions if alignment of individual subjects in the original photographs is not perfect. In fact, such distortions, known as “bending”, are common13,29,30,31 and apparent also in our data in both PC2 and PC3. Yet, once used for mapping, several QTL were detected for both of these PCs. Since it is unlikely that such QTL would have been detected if the actual variation in landmark locations had been obfuscated by bending effects, we believe that the landmark coordinates we used were still biologically informative. However, even if the validity of mapping results for variation in PC2 and PC3 coordinates were questioned, bending is unlikely to have influenced variation in all (ten) anatomical morphometric traits that were also used in this study.
In conclusion, we mapped several large-effect QTL, with a fairly high-level of accuracy, to different linkage groups in the nine-spined stickleback genome using a large panel of SNP-markers. These QTL influence variation in morphological traits of ecological and evolutionary significance, and the documented effects of allelic substitution in the QTL loci we detected align with the expectations based on the divergence in trait means among pond and marine populations in the wild. Hence, the results have taken us one step farther towards identifying genomic regions underlying differentiation among nine-spined stickleback populations living in contrasting environments. Comparisons to previous QTL studies conducted in this and related stickleback species show that the large-effect QTL detected in this study do not correspond to those detected in previous studies. This suggests that adaption to similar selection pressures may have been acquired using different genetic mechanisms.
Methods
Fish collection, crossing, and rearing
Grandparental sticklebacks (F0) from the Baltic Sea (Helsinki, Finland; 60°13′N, 25°11′E) and a pond (Rytilampi, Finland; 66°23′N, 29°19′E) population were collected in 2006 and transported to the aquarium facilities of the University of Helsinki. Details of crossing and fish rearing methods can be found from earlier studies112,113,118. In brief, a randomly selected female from the marine population and a male from the pond population were crossed artificially in July 2006. After F1 offspring from this cross reached sexual maturity, one randomly selected male and female were mated repeatedly to generate a F2 full-sib family comprised of seven successive clutches. F2 progeny were raised individually in 1.4 L tanks in zebrafish rack systems (Aquaneering Inc., San Diego, USA) at 17 °C for 187 days after hatching and then killed with an overdose of Tricaine Methanesulfonate (MS-222). Larval and juvenile sticklebacks were fed with brine shrimp nauplii and adults with frozen bloodworms. A total of 283 F2 progeny of seven clutches were available. To obtain morphological data, the F2 individuals were photographed with a digital camera for a lateral view. A fin clip from each individual was stored in 96% ethanol for DNA sequencing. The specimens were then fixed in 4% formalin for further phenotypic measurements. In addition, wild marine (16 individuals) and pond fish (63 individuals) used as reference material in this study were the same as used in Herczeg et al.35, but re-measured by the same person (JY) who also measured the F2 offspring to ensure comparable measurements. The experiments were conducted under the license from the Finnish National Animal Experiment Board (#STH379A). All experimental protocols were approved by the ethics committee of University of Helsinki, and all experiments were performed in accordance with relevant guidelines and regulations.
Morphological analysis
To quantify body shape variation, 21 landmarks (upper panel, Fig. 1) used in previous studies of three-spined13,29,119 and nine-spined sticklebacks35 were recorded from all F2 individuals using tpsDig v2.17120. MorphoJ121 was used to conduct geometric morphometric analyses of body shape variation. First, all the 21 landmarks were scaled and aligned through Procrustes superimpositions, and covariance matrix from the shape data was generated. Principal component analysis (PCA) was then applied to quantify shape variation among F2 progeny from the covariance matrix, working with partial warp scores. Since four degrees of freedom (cf. elimination of size, x and y coordinate positions, and orientation of specimens) were eliminated during Procrustes superimpositions122, variation along 38 (2k–4 = 38, where k = number of landmarks) PC axes were analyzed (see Results). PC scores were then used as separate traits to map QTL for shape variation. Although several earlier studies13,27,123,124 have used coordinates to map shape variation, our choice of using PC-scores was based on the fact that QTL inference based on shape coordinates can yield invalid inference125. Centroid size was also calculated with MorphoJ and used as covariate to control allometric effects during QTL mapping of PC scores. Apart from measurement of PC scores for the body shape, we also extracted information about variation in 11 anatomical measures (viz. head length, upper jaw length, lower jaw length, orbit diameter length, dorsal fin base length, anal fin base length, caudal peduncle length, caudal peduncle width, body depth, snout length, and standard body length) as defined in Fig. 1 (lower panel). In addition, the number of lateral plates on both sides of the body was counted from digital photographs. Each trait was measured or counted twice by the same person (JY), and the repeatabilities (R) of all measurements were high (R ≥ 0.89, P < 0.001 for all traits). Mean values of repeated measures were used in all analyses to minimize effect of measurement error. For QTL-mapping, the anatomical measurements were adjusted to variation in body size by regressing the trait values against standard body length and analyzing the residual variation. Sex of all F2 offspring was identified by gonadal inspection, and confirmed by genotyping a sex-linked microsatellite marker Stn19126.
Genotyping and linkage map construction
Genomic DNA was extracted from ethanol preserved fin clips using the phenol-chloroform method127. RAD library construction and sequencing were performed by BGI HONGKONG CO., LIMITED. Briefly, DNA was fragmented by the restriction enzyme PstI and DNA fragments of 300 to 500 bp were gel purified. Illumina sequencing adaptors and library specific barcodes were ligated to the digested DNA fragments and barcoded RAD samples were then pooled and sequenced on 24 lanes of the Illumina HiSeq2000 platform with 45 bp single end strategy. Adapters and barcodes were eliminated from reads and quality was checked using FastQC (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/).
Details of SNP calling and linkage map construction for this data have been described in Rastas et al.56. In brief, the reads of parental individuals without any missing nucleotides (Ns) were pooled and identical reads were grouped together, and only reads occurring between 10 and 1500 times were kept and taken as reference contigs; reads of each individual were mapped against the reference contigs with BWA128, together with SAMtools128, producing a single bam file for each individual. Individual genotype posterior probabilities, taking into account the read and mapping qualities, were obtained from the output of ‘mpileup’ on the bam files. Only positions with two or more alleles and at least three reads for more than 158 individuals were considered. The parental genotypes were called by maximizing the likelihood of offspring and parent genotypes and then each offspring was called with respect to its called parental genotypes. Parental genotypes were called only if their likelihood was 100 times higher than the second best parental genotype combination and the offspring genotypes were called similarly as the parents. We are aware of that genotyping with RAD-seq could lead to biases, such as those introduced by mutations in cut sites129. However, as we used a F2 recombinant cross and removed RAD-tags with extremely low or high read numbers (see above), such biases are unlikely to have had any large effect on our results.
Linkage group (LG) assignment was obtained using the software Lep-MAP256. First, SeparateChromosomes module was executed with LOD score limit 20 and minimum LG size of 10. Second, singular markers were added to the found LGs using JoinSingles module with LOD score limit of 10. Markers with more than 40 missing genotypes were removed from the LGs. Lep-MAP2 filters out markers by comparing the offspring genotype distribution and the expected Mendelian proportions with the dataTolerance of 0.01. Hence, all markers showing segregation distortion were filtered out and did not enter in the linkage map. The marker order was found using Lep-MAP2 allowing different recombination probabilities in both sexes. Ten independent runs were run and the marker order with best score56 was kept. Finally, markers with genotype error rate >0.1 were removed as well as markers from the map ends if they contributed over 10 cM (per marker) to map length. Likewise, if the parental coverage was above 500 (likely repetitive sequence) or below 20 (likely haplotype), these markers were removed.
QTL-mapping
Due to the lack of nine-spined stickleback genome sequence, we annotated the linkage map by utilizing the three-spined stickleback genome resource (Ensembl release-75). We used BLAST search with an E-value cutoff at 1 × 10−5 based on the genomic divergence between nine- and three-spined sticklebacks130 to map the nine-spined stickleback reference sequences with SNPs onto the three-spined stickleback genome. SNPs that could be uniquely mapped on the three-spined stickleback genome are defined as “informative SNPs”.
QTL-mapping analyses were performed with a two-stage approach (i.e. coarse+fine mapping) to overcome computational limitations in QTL mapping with a large number of markers131. To identify QTL across the genome, we initially conducted coarse-mapping using a simplified linkage map consisting of 466 randomly selected informative SNPs at intervals of approximately 5 cM. After the coarse-mapping, fine-mapping was performed with additional SNPs located near/in the detected QTL region(s) to provide better accuracy for QTL positions. The additional SNPs were selected based on a linkage map with 14,998 unique SNPs56 and added to the simplified linkage map to create a high-density linkage map for the fine-mapping.
QTL mapping was conducted using MapQTL 6.0131 based on interval mapping (IM) with the following parameters: mapping step size of five, maximum of 200 interactions, and a functional tolerance value of 1 × 10−8. Sex was included as a factor to control the gender dependent variation in the mapped traits. Significance of LOD thresholds were determined with 10 000 permutations in Permutation Test (PT) process using a significance level of P < 0.05. To confirm the IM results and to detect additional QTL, Multiple QTL-mapping (MQM) analyses were conducted using suggestive QTL as cofactors. Newly detected significant QTL were added as cofactors, and a new round of MQM mapping was then performed. This was repeated until no new significant QTL were detected128. Based on the population variance among F2 progeny, the proportion of the phenotypic variance explained (PVE) by a QTL was calculated in MapQTL 6.0. These procedures were followed both for coarse- and fine-mapping rounds. 1.5-LOD support intervals were calculated to provide 95% confidence interval for each QTL14,132.
SNP annotation and gene ontology analysis
The three-spined stickleback genome annotations were downloaded from Ensembl (release-75). BEDTools 2.17.0133 was used for annotation of the subset of SNPs identified by QTL-mapping to see whether the SNPs were located within or in proximity of known genes. Since the confidence intervals around detected QTL contained markers that we were unable to locate on the three-spined stickleback genome, we opted for the following approach. In brief, genomic intervals between the closest informative marker to the given QTL and the ends of its QTL confidence intervals were located in the three-spined stickleback genome. All the genes within these genomic intervals were then listed and candidate genes for given QTL were identified on the basis of gene ontology information. If gene ontology information did not retrieve any obvious candidate genes, none were called (but the identities of the genes within the search interval are reported). Gene ontology (GO) terms of the three-spined stickleback genes were retrieved with BioMart134 from Ensembl.
Data accessibility
All data used in this study was given in the supplementary materials. The linkage map and genotype data are given in Supplementary Table 1. The original measurements of morphological traits in the F2 progeny are given in Supplementary Table 2.
Additional Information
How to cite this article: Yang, J. et al. Quantitative trait locus analysis of body shape divergence in nine-spined sticklebacks based on high-density SNP-panel. Sci. Rep. 6, 26632; doi: 10.1038/srep26632 (2016).
References
Bickford, D. et al. Cryptic species as a window on diversity and conservation. Trends Ecol. Evol. 22, 148–155 (2007).
Fruciano, C., Tigano, C. & Ferrito, V. Geographical and morphological variation within and between colour phases in Coris julis (L. 1758), a protogynous marine fish. Biol. J. Linn. Soc. 104, 148–162 (2011).
Foster, K., Bower, L. & Piller, K. Getting in shape: habitat-based morphological divergence for two sympatric fishes. Biol. J. Linn. Soc. 114, 152–162 (2015).
Hendry, A. P., Kelly, M. L., Kinnison, M. T. & Reznick, D. N. Parallel evolution of the sexes? Effects of predation and habitat features on the size and shape of wild guppies. J. Evol. Biol. 19, 741–754 (2006).
Schluter, D. Adaptive radiation along genetic lines of least resistance. Evolution 50, 1766–1774 (1996).
Leinonen, T., McCairns, R. J. S., O’Hara, R. B. & Merilä, J. Q ST –F ST comparisons: evolutionary and ecological insights from genomic heterogeneity. Nat. Rev. Genet. 14, 179–190 (2013).
Karhunen, M., Ovaskainen, O., Herczeg, G. & Merilä, J. Bringing habitat information into statistical tests of local adaptation in quantitative traits: a case study of nine-spined sticklebacks. Evolution 68, 559–568 (2014).
McCarthy, M. I. et al. Genome-wide association studies for complex traits: consensus, uncertainty and challenges. Nat. Rev. Genet. 9, 356–369 (2008).
Solovieff, N., Cotsapas, C., Lee, P. H., Purcell, S. M. & Smoller, J. W. Pleiotropy in complex traits: challenges and strategies. Nat. Rev. Genet. 14, 483–495 (2013).
Burns, J. G., Di Nardo, P. & Rodd, F. H. The role of predation in variation in body shape in guppies Poecilia reticulata: a comparison of field and common garden phenotypes. J. Fish Biol. 75, 1144–1157 (2009).
Leinonen, T., Cano, J. M., Mäkinen, H. & Merilä, J. Contrasting patterns of body shape and neutral genetic divergence in marine and lake populations of threespine sticklebacks. J. Evol. Biol. 19, 1803–1812 (2006).
Hassell, E. M., Meyers, P. J., Billman, E. J., Rasmussen, J. E. & Belk, M. C. Ontogeny and sex alter the effect of predation on body shape in a livebearing fish: sexual dimorphism, parallelism, and costs of reproduction. Ecol. Evol. 2, 1738–1746 (2012).
Albert, A. Y. et al. The genetics of adaptive shape shift in stickleback: pleiotropy and effect size. Evolution 62, 76–85 (2008).
Franchini, P. et al. Genomic architecture of ecologically divergent body shape in a pair of sympatric crater lake cichlid fishes. Mol. Ecol. 23, 1828–1845 (2014).
Mackay, T. F., Stone, E. A. & Ayroles, J. F. The genetics of quantitative traits: challenges and prospects. Nat. Rev. Genet. 10, 565–577 (2009).
Lynch, M. & Walsh, B. Genetics and analysis of quantitative traits, 1–316 (Sinauer Associates, Sunderland, 1998).
Reid, D. T. & Peichel, C. L. Perspectives on the genetic architecture of divergence in body shape in sticklebacks. Integr. Comp. Biol. 50, 1057–1066 (2010).
Bell, M. A., Stewart, J. D. & Park, P. J. The world’s oldest fossil threespine stickleback fish. Copeia 2009, 256–265 (2009).
Merilä, J. Nine-spined stickleback (Pungitius pungitius): an emerging model for evolutionary biology research. Ann. NY Acad. Sci. 1289, 18–35 (2013).
Wootton, R. J. The biology of the sticklebacks, 1–387 (Academic Press, London, 1976).
Wootton, R. J. A functional biology of sticklebacks, 1–265 (University of California Press, Berkeley, 1984).
Shapiro, M. D. et al. Genetic and developmental basis of evolutionary pelvic reduction in threespine sticklebacks. Nature 428, 717–723 (2004).
Peichel, C. L. et al. The genetic architecture of divergence between threespine stickleback species. Nature 414, 901–905 (2001).
Colosimo, P. F. et al. Widespread parallel evolution in sticklebacks by repeated fixation of ectodysplasin alleles. Science 307, 1928–1933 (2005).
Colosimo, P. F. et al. The genetic architecture of parallel armor plate reduction in threespine sticklebacks. PloS Biol. 2, e109 (2004).
Schluter, D., Marchinko, K. B., Barrett, R. D. & Rogers, S. M. Natural selection and the genetics of adaptation in threespine stickleback. Philos. Trans. R. Soc. B 365, 2479–2486 (2010).
Liu, J. et al. Identification of major and minor QTL for ecologically important morphological traits in three-spined sticklebacks (Gasterosteus aculeatus). G3-Genes Genom. Genet. 4, 595–604 (2014).
Malek, T. B., Boughman, J. W., Dworkin, I. & Peichel, C. L. Admixture mapping of male nuptial colour and body shape in a recently formed hybrid population of threespine stickleback. Mol. Ecol. 21, 5265–5279 (2012).
Rogers, S. M. et al. Genetic signature of adaptive peak shift in threespine stickleback. Evolution 66, 2439–2450 (2012).
Leinonen, T., Cano, J. M. & Merilä, J. Genetics of body shape and armour variation in threespine sticklebacks. J. Evol. Biol. 24, 206–218 (2011).
Herczeg, G., Turtiainen, M. & Merilä, J. Morphological divergence of North-European nine-spined sticklebacks (Pungitius pungitius): signatures of parallel evolution. Biol. J. Linn. Soc. 101, 403–416 (2010).
Gonda, A., Herczeg, G. & Merilä, J. Adaptive brain size divergence in nine-spined sticklebacks (Pungitius pungitius)? J. Evol. Biol. 22, 1721–1726 (2009).
Herczeg, G., Gonda, A. & Merilä, J. Predation mediated population divergence in complex behaviour of nine-spined stickleback (Pungitius pungitius). J. Evol. Biol. 22, 544–552 (2009).
Herczeg, G., Gonda, A. & Merilä, J. Evolution of gigantism in nine-spined sticklebacks. Evolution 63, 3190–3200 (2009).
Shimada, Y. et al. Quantitative genetics of body size and timing of maturation in two nine-spined stickleback (Pungitius pungitius) populations. PloS One 6, e28859 (2011).
Ab Ghani, N. I., Herczeg, G. & Merilä, J. Body size divergence in nine-spined sticklebacks: disentangling additive genetic and maternal effects. Biol. J. Linn. Soc. 107, 521–528 (2012).
Välimäki, K., Herczeg, G. & Merilä, J. Morphological anti-predator defences in the nine-spined stickleback: constitutive, induced or both? Biol. J. Linn. Soc. 107, 854–866 (2012).
Bradshaw, W. E., Emerson, K. J., Catchen, J. M., Cresko, W. A. & Holzapfel, C. M. Footprints in time: comparative quantitative trait loci mapping of the pitcher-plant mosquito, Wyeomyia smithii . Proc. R. Soc. B 279, 4551–4558 (2012).
Bernatchez, L. et al. On the origin of species: insights from the ecological genomics of lake whitefish. Philos. Trans. R. Soc. B 365, 1783–1800 (2010).
Stapley, J. et al. Adaptation genomics: the next generation. Trends Ecol. Evol. 25, 705–712 (2010).
Mardis, E. R. Next-generation DNA sequencing methods. Annu. Rev. Genomics Hum. Genet. 9, 387–402 (2008).
Imelfort, M., Duran, C., Batley, J. & Edwards, D. Discovering genetic polymorphisms in next-generation sequencing data. Plant Biotechnol. J. 7, 312–317 (2009).
Luikart, G., England, P. R., Tallmon, D., Jordan, S. & Taberlet, P. The power and promise of population genomics: From genotyping to genome typing. Nat. Rev. Genet. 4, 981–994 (2003).
Andolfatto, P. et al. Multiplexed shotgun genotyping for rapid and efficient genetic mapping. Genome Res. 21, 610–617 (2011).
Barbazuk, W. B., Bedell, J. A. & Rabinowicz, P. D. Reduced representation sequencing: a success in maize and a promise for other plant genomes. Bioessays 27, 839–848 (2005).
Van Tassell, C. P. et al. SNP discovery and allele frequency estimation by deep sequencing of reduced representation libraries. Nat. Methods 5, 247–252 (2008).
Baird, N. A. et al. Rapid SNP discovery and genetic mapping using sequenced RAD markers. PloS One 3, e3376 (2008).
Davey, J. W. et al. Genome-wide genetic marker discovery and genotyping using next-generation sequencing. Nat. Rev. Genet. 12, 499–510 (2011).
Nielsen, R., Paul, J. S., Albrechtsen, A. & Song, Y. S. Genotype and SNP calling from next-generation sequencing data. Nat. Rev. Genet. 12, 443–451 (2011).
Baxter, S. W. et al. Linkage mapping and comparative genomics using next-generation RAD sequencing of a non-model organism. PloS One 6, e19315 (2011).
Chutimanitsakun, Y. et al. Construction and application for QTL analysis of a Restriction Site Associated DNA (RAD) linkage map in barley. BMC Genomics 12, 4 (2011).
Pfender, W. F., Saha, M. C., Johnson, E. A. & Slabaugh, M. B. Mapping with RAD (restriction-site associated DNA) markers to rapidly identify QTL for stem rust resistance in Lolium perenne . Theor. Appl. Genet. 122, 1467–1480 (2011).
Würschum, T. & Mapping QTL for agronomic traits in breeding populations. Theor. Appl. Genet. 125, 201–210 (2012).
Stange, M., Utz, H. F., Schrag, T. A., Melchinger, A. E. & Würschum, T. High-density genotyping: an overkill for QTL mapping? Lessons learned from a case study in maize and simulations. Theor. Appl. Genet. 126, 2563–2574 (2013).
Ocalewicz, K., Fopp-Bayat, D., Woznicki, P. & Jankun, M. Heteromorphic sex chromosomes in the ninespine stickleback Pungitius pungitius . J. Fish Biol. 73, 456–462 (2008).
Rastas, P., Calboli, F. C. F., Guo, B., Shikano, T. & Merilä, J. Construction of ultra-dense linkage maps with Lep-MAP2: stickleback F2 recombinant crosses as an example. Genome Biol. Evol. 8, 78–93 (2016).
Cano, J. M., Matsuba, C., Makinen, H. & Merilä, J. The utility of QTL-Linked markers to detect selective sweeps in natural populations-a case study of the EDA gene and a linked marker in threespine stickleback. Mol. Ecol. 15, 4613–4621 (2006).
Shapiro, M. D. et al. The genetic architecture of skeletal convergence and sex determination in ninespine sticklebacks. Curr. Biol. 19, 1156–1156 (2009).
Keivany, Y. & Nelson, J. S. Taxonomic review of the genus Pungitius, ninespine sticklebacks (Gasterosteidae). Cybium 24, 107–122 (2000).
Mattern, M. Y. In Biology of the Three-spined Sticklebacks (eds Östlund-Nilsson, S., Mayer, I. & Huntingford, F. A. ) 1–40 (CRC Press, Boca Raton, 2007).
Videler, J. J. On the interrelationships between morphology and movement in the tail of the cichlid fish Tilapia nilotica (L.). Neth. J. Zool. 25, 143–194 (1974).
Gibb, A. C., Jayne, B. C. & Lauder, G. V. Kinematics of pectoral fin locomotion in the bluegill sunfish Lepomis macrochirus . J. Exp. Biol. 189, 133–161 (1994).
Fish, F. E. Balancing requirements for stability and maneuverability in cetaceans. Integr. Comp. Biol. 42, 85–93 (2002).
Blake, R. W. Fish functional design and swimming performance. J. Fish Biol. 65, 1193–1222 (2004).
Videler, J. J. Fish swimming, 73–74 (Chapman and Hall, London, 1993).
Leinonen, T., Herczeg, G., Cano, J. M. & Merilä, J. Predation-imposed selection on threespine stickleback (Gasterosteus aculeatus) morphology: a test of the refuge use hypothesis. Evolution 65, 2916–2926 (2011).
Barrett, R. D. H. Adaptive evolution of lateral plates in three-spined stickleback Gasterosteus aculeatus: a case study in functional analysis of natural variation. J. Fish Biol. 77, 311–328 (2010).
Barrett, R. D. H., Rogers, S. M. & Schluter, D. Natural selection on a major armor gene in threespine stickleback. Science 322, 255–257 (2008).
Kitano, J. et al. Reverse evolution of armor plates in the threespine stickleback. Curr. Biol. 18, 769–774 (2008).
Bell, M. A. Lateral plate evolution in the threespine stickleback: getting nowhere fast. Genetica 112, 445–461 (2001).
Cresko, W. A. et al. Parallel genetic basis for repeated evolution of armor loss in Alaskan threespine stickleback populations. Proc. Natl. Acad. Sci. USA 101, 6050–6055 (2004).
DeFaveri, J. & Merilä, J. Evidence for adaptive phenotypic differentiation in Baltic Sea sticklebacks. J. Evol. Biol. 26, 1700–1715 (2013).
Dowling, T. E., Martasian, D. P. & Jeffery, W. R. Evidence for multiple genetic forms with similar eyeless phenotypes in the blind cavefish, Astyanax mexicanus. Mol. Biol. Evol. 19, 446–455 (2002).
Merilä, J. Lakes and ponds as model systems to study parallel evolution. J. Limnol. 73, 33–45 (2014).
Foolad, M. R. & Chen, F. Q. RFLP mapping of QTLs conferring salt tolerance during the vegetative stage in tomato. Theor. Appl. Genet. 99, 235–243 (1999).
Collard, B. C. Y., Jahufer, M. Z. Z., Brouwer, J. B. & Pang, E. C. K. An introduction to markers, quantitative trait loci (QTL) mapping and marker-assisted selection for crop improvement: The basic concepts. Euphytica 142, 169–196 (2005).
Barchi, L. et al. A RAD tag derived marker based eggplant linkage map and the location of QTLs determining anthocyanin pigmentation. PloS One 7, e43740 (2012).
Slate, J. et al. A genome scan for quantitative trait loci in a wild population of red deer (Cervus elaphus). Genetics 162, 1863–1873 (2002).
Slate, J. From beavis to beak color: a simulation study to examine how much QTL mapping can reveal about the genetic architecture of quantitative traits. Evolution 67, 1251–1262 (2013).
Beavis, W. D. In Proceedings of the Forty-Ninth Annual Corn & Sorghum industry research conference (ed Wilkinson, D. B. ) 250–266 (American Seed Trade Association, Washington D. C., 1994).
Beavis, W. D. In Molecular Dissection of Complex Traits (ed. Paterson, A. H. ) 145–162 (CRC Press, Boca Raton, 1998).
Pallares, L. F., Harr, B., Turner, L. M. & Tautz, D. Use of a natural hybrid zone for genomewide association mapping of craniofacial traits in the house mouse. Mol. Ecol. 23, 5756–5770 (2014).
Gamazon, E. R., Cox, N. J. & Davis, L. K. Structural architecture of SNP effects on complex traits. Am. J. Hum. Genet. 95, 477–489 (2014).
Polderman, T. J. C. et al. Meta-analysis of the heritability of human traits based on fifty years of twin studies. Nat. Genet. 47, 702–709 (2015).
Bell, M. A., Orti, G., Walker, J. A. & Koenings, J. P. Evolution of pelvic reduction in threespine stickleback fish-a test of competing hypotheses. Evolution 47, 906–914 (1993).
DeFaveri, J., Viitaniemi, H., Leder, E. & Merilä, J. Characterizing genic and nongenic molecular markers: comparison of microsatellites and SNPs. Mol. Ecol. Resour. 13, 377–392 (2013).
Shimada, Y., Shikano, T. & Merilä, J. A high incidence of selection on physiologically important genes in the three-spined stickleback, Gasterosteus aculeatus. Mol. Biol. Evol. 28, 181–193 (2011).
Galeano, C. H. et al. Gene-based single nucleotide polymorphism markers for genetic and association mapping in common bean. BMC Genet. 13, 48 (2012).
Vilas, R., Bouza, C., Vera, M., Millán, A. & Martínez, P. Variation in anonymous and EST-microsatellites suggests adaptive population divergence in turbot. Mar. Ecol. Prog. Ser. 420, 231–239 (2010).
Cortés, A. J., Chavarro, M. C. & Blair, M. W. SNP marker diversity in common bean (Phaseolus vulgaris L.). Theor. Appl. Genet. 123, 827–845 (2011).
Sabherwal, N. et al. Long-range conserved non-coding SHOX sequences regulate expression in developing chicken limb and are associated with short stature phenotypes in human patients. Hum. Mol. Genet. 16, 210–222 (2007).
Klopocki, E. & Mundlos, S. Copy-number variations, noncoding sequences, and human phenotypes. Annu. Rev. Genomics Hum. Genet. 12, 53–72 (2011).
Vockley, C. M. et al. Massively parallel quantification of the regulatory effects of noncoding genetic variation in a human cohort. Genome Res. 25, 1206–1214 (2015).
Melvin, V. S., Feng, W., Hernandez-Lagunas, L., Artinger, K. B. & Williams, T. A morpholino-based screen to identify novel genes involved in craniofacial morphogenesis. Dev. Dynam. 242, 817–831 (2013).
Powder, K. E., Cousin, H., McLinden, G. P. & Albertson, R. C. A nonsynonymous mutation in the transcriptional regulator lbh is associated with cichlid craniofacial adaption and neural crest cell development. Mol. Biol. Evol. 31, 3113–3124 (2014).
Hall, C., Flores, M. V., Murison, G., Crosier, K. & Crosier, P. An essential role for zebrafish Fgfrl1 during gill cartilage development. Mech. Develop. 123, 925–940 (2006).
Gessert, S., Bugner, V., Tecza, A., Pinker, M. & Kuhl, M. FMR1/FXR1 and the miRNA pathway are required for eye and neural crest development. Dev. Biol. 341, 222–235 (2010).
Murciano, C. et al. Ray-interray interactions during fin regeneration of Danio rerio . Dev. Biol. 252, 214–224 (2002).
Hammerschmidt, M., Serbedzija, G. N. & McMahon, A. P. Genetic analysis of dorsoventral pattern formation in the zebrafish: requirement of a BMP-like ventralizing activity and its dorsal repressor. Genes Dev. 10, 2452–2461 (1996).
Vauti, F. et al. The mouse Trm1-like gene is expressed in neural tissues and plays a role in motor coordination and exploratory behaviour. Gene 389, 174–185 (2007).
Galligan, J. J. & Petersen, D. R. The human protein disulfide isomerase gene family. Hum. Genomics 6, 6 (2012).
Laine, V. N., Shikano, T., Herczeg, G., Vilkki, J. & Merilä, J. Quantitative trait loci for growth and body size in the nine-spined stickleback Pungitius pungitius L. Mol. Ecol. 22, 5861–5876 (2013).
Laine, V. N., Herczeg, G., Shikano, T., Vilkki, J. & Merilä, J. QTL analysis of behavior in nine-spined sticklebacks (Pungitius pungitius). Behav. Genet. 44, 77–88 (2014).
Makinen, H. S., Cano, J. M. & Merilä, J. Identifying footprints of directional and balancing selection in marine and freshwater three-spined stickleback (Gasterosteus aculeatus) populations. Mol. Ecol. 17, 3565–3582 (2008).
Miller, C. T. et al. Modular skeletal evolution in sticklebacks is controlled by additive and clustered quantitative trait loci. Genetics 197, 405–420 (2014).
Greenwood, A. K., Wark, A. R., Yoshida, K. & Peichel, C. L. Genetic and neural modularity underlie the evolution of schooling behavior in threespine sticklebacks. Curr. Biol. 23, 1884–1888 (2013).
Greenwood, A. K. et al. Genetic mapping of natural variation in schooling tendency in the threespine stickleback. G3-Genes Genom. Genet. 5, 761–769 (2015).
Wu, K. et al. High-density genetic map construction and QTLs analysis of grain yield-related traits in Sesame (Sesamum indicum L.) based on RAD-Seq techonology. BMC Plant Biol. 14, 274 (2014).
Li, Y. G. & He, M. X. Genetic mapping and QTL analysis of growth-related traits in Pinctada fucata using restriction-site associated DNA sequencing. PloS One 9, e111707 (2014).
Houston, R. D. et al. Characterisation of QTL-linked and genome-wide restriction site-associated DNA (RAD) markers in farmed Atlantic salmon. BMC Genomics 13, 244 (2012).
Hecht, B. C., Thrower, F. P., Hale, M. C., Miller, M. R. & Nichols, K. M. Genetic architecture of migration-related traits in rainbow and steelhead trout, Oncorhynchus mykiss . G3-Genes Genom. Genet. 2, 1113–1127 (2012).
Palaiokostas, C. et al. Mapping the sex determination locus in the Atlantic halibut (Hippoglossus hippoglossus) using RAD sequencing. BMC Genomics 14, 566 (2013).
Gagnaire, P. A., Normandeau, E., Pavey, S. A. & Bernatchez, L. Mapping phenotypic, expression and transmission ratio distortion QTL using RAD markers in the lake whitefish (Coregonus clupeaformis). Mol. Ecol. 22, 3036–3048 (2013).
Kakioka, R., Kokita, T., Kumada, H., Watanabe, K. & Okuda, N. A RAD-based linkage map and comparative genomics in the gudgeons (genus Gnathopogon, Cyprinidae). BMC Genomics 14, 32 (2013).
Gonen, S. et al. Linkage maps of the Atlantic salmon (Salmo salar) genome derived from RAD sequencing. BMC Genomics 15, 166 (2014).
Berner, D., Moser, D., Roesti, M., Buescher, H. & Salzburger, W. Genetic architecture of skeletal evolution in European lake and stream stickleback. Evolution 68, 1792–1805 (2014).
Mitteroecker, P. & Bookstein, F. Linear discrimination, ordination, and the visualization of selection gradients in modern morphometrics. Evol. Biol. 38, 100–114 (2011).
Shikano, T., Laine, V. N., Herczeg, G., Vilkki, J. & Merilä, J. Genetic architecture of parallel pelvic reduction in ninespine sticklebacks. G3-Genes Genom. Genet. 3, 1833–1842 (2013).
Schluter, D., Clifford, E. A., Nemethy, M. & McKinnon, J. S. Parallel evolution and inheritance of quantitative traits. Am. Nat. 163, 809–822 (2004).
Rohlf, F. J. TpsDig v. 2.17. Department of Ecology and Evolution, State University of New York at Stony Brook, Stony Brook, New York, USA URL http://life.bio.sunysb.edu/morph/ (2013).
Klingenberg, C. P. MorphoJ: an integrated software package for geometric morphometrics. Mol. Ecol. Resour. 11, 353–357 (2011).
Klingenberg, C. P., Leamy, L. J., Routman, E. J. & Cheverud, J. M. Genetic architecture of mandible shape in mice: effects of quantitative trait loci analyzed by geometric morphometrics. Genetics 157, 785–802 (2001).
Laporte, M. et al. RAD-QTL mapping reveals both genome-level parallelism and different genetic architecture underlying the evolution of body shape in lake whitefish (Coregonus clupeaformis) species pairs. G3-Genes Genom. Genet. 5, 1481–1491 (2015).
Kakioka, R., Kokita, T., Kumada, H., Watanabe, K. & Okuda, N. Genomic architecture of habitat-related divergence and signature of directional selection in the body shapes of Gnathopogon fishes. Mol. Ecol. 24, 4159–4174 (2015).
Márquez, E. J., Cabeen, R., Woods, R. P. & Houle, D. The measurement of local variation in shape. Evol. Biol. 39, 419–439 (2012).
Shikano, T., Herczeg, G. & Merilä, J. Molecular sexing and population genetic inference using a sex-linked microsatellite marker in the nine-spined stickleback (Pungitius pungitius). BMC Res . Notes 4, 119 (2011).
Taggart, J. B., Hynes, R. A., Prodöhl, P. A. & Ferguson, A. A simplified protocol for routine total DNA isolation from salmonid fishes. J. Fish Biol. 40, 963–965 (1992).
Li, H. et al. The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079 (2009).
Gautier, M. et al. The effect of RAD allele dropout on the estimation of genetic variation within and between populations. Mol. Ecol. 22, 3165–3178 (2013).
Guo, B. C., Chain, F. J. J., Bornberg-Bauer, E., Leder, E. H. & Merilä, J. Genomic divergence between nine- and three-spined sticklebacks. BMC Genomics 14, 756 (2013).
Van Ooijen, J. W. MapQTL v. 6.0: Software for the mapping of quantitative trait loci in experimental populations of diploid species. Kyazma BV, Wageningen, The Netherlands. UTL https://www.kyazma.nl/index.php/MapQTL/ (2009).
Dupuis, J. & Siegmund, D. Statistical methods for mapping quantitative trait loci from a dense set of markers. Genetics 151, 373–386 (1999).
Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841–842 (2010).
Smedley, D. et al. BioMart-biological queries made easy. BMC Genomics 10, 22 (2009).
Acknowledgements
We thank Abigél Gonda, Gabor Herczeg and Mirva Turtiainen for help with the fish breeding, Kirsi Kähkönen and Marika Karjalainen for assistance in the laboratory, Johan W. van Ooijen for technical assistance in QTL-mapping analyses, and Federico Calboli and Jacquelin De Faveri for comments on earlier version of this manuscript. We also thank the computing resource support from CSC (the Finnish IT Center for Science Ltd administered by the Ministry of Education and Culture, Finland). Our research was supported by the Academy of Finland (grants 129662, 134728 and 218343 to J.M.) and by the China Scholarship Council (J.Y.).
Author information
Authors and Affiliations
Contributions
J.M. conceived the study. T.S. and J.M. carried out the laboratory experiments. J.Y. and B.G. analyzed the data. J.Y., B.G. and J.M. wrote the manuscript with contributions from T.S. and X.L. All authors reviewed the manuscript.
Corresponding authors
Ethics declarations
Competing interests
The authors declare no competing financial interests.
Supplementary information
Rights and permissions
This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/
About this article
Cite this article
Yang, J., Guo, B., Shikano, T. et al. Quantitative trait locus analysis of body shape divergence in nine-spined sticklebacks based on high-density SNP-panel. Sci Rep 6, 26632 (2016). https://doi.org/10.1038/srep26632
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/srep26632
This article is cited by
-
Quantitative trait loci for morphometric traits in multiple families of common carp (Cyprinus carpio)
Science China Life Sciences (2017)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.
