Causes and consequences of intra-specific variation in vertebral number

Intraspecific variation in vertebral number is taxonomically widespread. Much scientific attention has been directed towards understanding patterns of variation in vertebral number among individuals and between populations, particularly across large spatial scales and in structured environments. However, the relative role of genes, plasticity, selection, and drift as drivers of individual variation and population differentiation remains unknown for most systems. Here, we report on patterns, causes and consequences of variation in vertebral number among and within sympatric subpopulations of pike (Esox lucius). Vertebral number differed among subpopulations, and common garden experiments indicated that this reflected genetic differences. A QST-FST comparison suggested that population differences represented local adaptations driven by divergent selection. Associations with fitness traits further indicated that vertebral counts were influenced both by stabilizing and directional selection within populations. Overall, our study enhances the understanding of adaptive variation, which is critical for the maintenance of intraspecific diversity and species conservation.

Scientific RepoRts | 6:26372 | DOI: 10.1038/srep26372 and juveniles (χ 2 = 27.57, df = 6, N = 146, P < 0.0001, Fig. 2). Despite that breeding streams were only separated by short distances (Fig. 1a), the magnitude of the differences between population averages was large, corresponding to 18% (1.4/8, adults) and 14% (1.0/7, juveniles) of the total range of variation seen among individuals, and to 1. Variation in vertebral number among subpopulations accounted for by genetic modifications. A central question was whether differences among subpopulations in VN reflected genetic modifications, developmental plasticity or a combination of the two. This was evaluated using data for juveniles reared in the common garden experiment, the results of which uncovered significant differences in frequency distribution of VN among subpopulations (Fisher's Exact test: χ 2 = 56.41, df = 6, P < 0.0001, Fig. 2c). The magnitude of the differences between subpopulation averages in the common garden corresponded to 18% (0.9/5) of the total range of variation seen among individuals, and to 1.1 (0.9/0.8) of the standard deviation of the pooled sample of individuals (Fig. 2). This result indicated that differences were, at least in part, genetically determined.
Resemblance between captive reared and wild-caught juveniles indicated that population differences reflected genetic modifications rather than plasticity. In two subpopulations (L & O) data allowed for comparisons of VN between wild-caught and captive-reared (common garden) juveniles to evaluate the hypothesis that differences between subpopulations represented plastic responses to environmental components. However, VN did not differ between juveniles that had developed in different environments in any of the subpopulations (subpopulation L: Fisher's Exact test: χ 2 = 3.97, df = 3, P = 0.22; subpopulation O: Fisher's Exact test: χ 2 = 3.15, df = 3, P = 0.38, Fig. 2). The range of VN in captive-reared juveniles (59)(60)(61)(62)(63) as well as the magnitude of the differences between subpopulations resembled those seen in wild-caught juveniles (Fig. 2). This similarity between wild-caught and common garden reared juveniles indicated that the variation among subpopulations seen in wild-caught individuals likely represented genetic modifications rather than developmental plasticity.
Developmental plasticity of vertebral number in response to incubation temperature. The split-brood experiment, in which full-siblings were exposed to different incubation temperatures (5,10, 15 and 20 °C), uncovered that VN decreased significantly with increasing temperatures (GLM, effect of treatment: F 2,36 = 4.82, P = 0.014, effect size η 2 = 0.09) and that there were significant effects of family (GLM, effect of family: F 4,36 = 5.89, P = 0.001, η 2 = 0.22) (Fig. 3). None of the eggs incubated in the 20 °C treatment hatched. The results from the split-brood experiment further uncovered crossing norms of reaction (as evidenced by a significant effect on vertebral number of the family by temperature interaction, GLM: F 7,36 = 3.01, P = 0.014, η 2 = 0.20, Fig. 3), suggesting that there existed genetic variation for plasticity.  Table 1.
Scientific RepoRts | 6:26372 | DOI: 10.1038/srep26372 Q ST -F ST comparison revealed that the divergence in vertebral number among subpopulations was driven by selection. The quantitative genetic differentiation (as estimated by Q ST ) of VN was high (Q ST = 0.54, 95% CI = 0.24-1.00) and it significantly exceeded the neutral genetic expectation set by F ST (0.069, 95% CI = 0.048-0.085), as illustrated by non-overlapping 95% confidence intervals. Theory states that subpopulation divergence is due to divergent selection if Q ST is significantly different from F  . Our results thus indicated that differences among subpopulations could be ascribed to divergent selection rather than to genetic drift.

Differences in vertebral number between adults and juveniles indicated viability selection.
Comparisons of VN frequency distributions revealed differences between wild-caught adults and juveniles, after controlling for the variation among subpopulations (Log-linear model; effect of life-stage: χ 2 = 17.85, df = 1, P < 0.0001; effect of subpopulation: χ 2 = 11.47, df = 2, P = 0.003). Adults had lower VN than juveniles in all three subpopulations (Fig. 2), suggesting that variability among individuals within populations was affected by natural (viability) selection.
Association of vertebral number with growth rate in juveniles and subadult life-stages. Data on body size indicated stabilizing selection 48 on VN in juveniles. Individuals with an intermediate number of vertebrae were larger on average compared with individuals having either relatively few or relatively many vertebrae (Fig. 4a,b). Quadratic regression demonstrated a statistically significant curvilinear relationship between body size and VN in wild-caught juveniles (y = − 4183 + 139.1 (± 57.6 SE) X-1.14 (± 0.475 SE) X 2 ; test of linear effect: t = 2.42, P = 0.0184, effect size η 2 = 0.08; test of quadratic effect: t = − 2.41, P = 0.0186, η 2 = 0.08, randomization test for significant quadratic effect, p = 0.029 (28/2000)) but not in common garden reared juveniles  (y = − 2757 + 92.1 (± 60.93) X-0.75 (± 0.50) X 2 ; test of linear effect: t = 1.51, P = 0.136, effect size η 2 = 0.03; test of quadratic effect: t = − 1.50, P = 0.138, η 2 = 0.03, randomization test for significant quadratic effect, p = 0.26 (259/2000)). These results indicated that the relationship linking body length to VN was expressed more strongly among individuals under natural conditions in the wild than in the laboratory (Fig. 4a,b). However, the shape of the relationship did not differ between captive-reared and wild-caught juveniles (effect of interaction between origin and the linear effect:

Discussion
We examined causes and consequences of variation in VN among and within six sympatric anadromous subpopulations of E. lucius pike using data for 721 wild-caught and captive-reared individuals. Our main findings were: i) number of vertebrae varied among subpopulations in the natural environment, ii) variation among subpopulations was, at least in part, genetically based and driven by divergent stage-specific selection as evidenced by a common garden experiment and a Q ST -F ST comparison, iii) resemblance between captive-reared and wild-caught individuals suggested that the pattern of variation documented in natural subpopulations reflected genetic differences rather than plasticity, iv) cross-sectional comparisons indicated directional viability selection on VN, v) associations with juvenile body size, growth trajectories and reproductive effort suggested that individuals with an intermediate number of vertebrae had higher fitness, indicative of stabilizing selection within subpopulations.
Population divergence in VN has typically been evaluated across large geographic distances, commonly to assess whether patterns of variation conform to Jordan's rule 17,19 . Considerably less attention has been directed towards documenting and understanding patterns of variation among populations within less delineated systems at small spatial scales 2,15,21 . In this study, we demonstrated that the frequency distribution of VN varied among pike that were sympatric in the Baltic Sea but originated from different subpopulations that utilize specific breeding streams due to homing behaviour 39 . The magnitude of the differences between subpopulation averages documented in this study (which corresponded to 1.26 s.d.) was large compared to previous, albeit few, studies of divergence in VN in confined systems and sympatric populations. For example, Manier, et al. 15 found that divergence among six local populations of the Garter snake (Thamnophis elegans) corresponded to 0.5 s.d., whereas Aguirre, et al. 9 found that two partially sympatric populations (anadromous and resident) of three-spine stickleback (Gasterosteus aculetus) were differentiated by 0.09 s.d.
Phenotypic variation can represent underlying genetic differences, environmentally induced plasticity, or reflect a combination of the two [50][51][52] . Several studies have shown that genetic components are involved in intraspecific variation of VN and that heritability is relatively high 2,4,8,28 . The results from our common garden experiment indicated that variation in VN among subpopulations was, at least in part, genetically based. Because we used only a F1-generation non-genetic maternal effects could also have been involved 4,8 . However, we produced common garden juveniles by stripping gametes from adults living in a shared environment prior to reproduction 39 . The influence of any environmentally induced (maternal effects) variation should therefore have been negligible.
Few studies have evaluated the relative roles of genetic components and plasticity in shaping variation in VN among populations 2,7,8 . It is well established that incubation temperature during embryogenesis can induce variation in VN 18,24 . To evaluate the role of temperature-induced plasticity we manipulated the incubation temperature in a split-brood experiment. Admittedly, the number of families in the split-brood experiment was small. Despite this, results uncovered that VN decreased with increasing temperature, which is coherent with most prevailing empirical evidence in other species 4,24 . Moreover, this experiment revealed genotype by environment effects on VN, emphasizing that there existed genetic variation of reaction norms allowing for evolution of developmental plasticity 53 . While it is possible that temperature-induced plasticity could have contributed to the differences in VN among subpopulations documented in wild-caught individuals, this would have required considerable and unrealistic temperature differences (approximately ∼ 10 °C, compare Figs 2 and 3) between breeding habitats. Any such temperature differences between breeding habitats should have been manifest as a disparity in VN between wild-caught and captive-reared juveniles. However, results for wild-caught and captive-reared juveniles were similar, suggesting that subpopulation differences reflected genetic effects rather than plastic responses to different breeding environments.
Population divergence in VN has been described in other species but less explored is whether such divergence is adaptive 15,21 . We addressed this by comparing phenotypic divergence (Q ST ) in the common garden experiment with divergence at neutral microsatellite loci (F ST ) and found that evolutionary divergence among subpopulations could be better explained by divergent selection than by random genetic drift. The Q ST -F ST comparison included only three of the six populations. The results nevertheless indicated that differences in VN among these three subpopulations were shaped by local adaptations to their specific breeding streams. That subpopulations were subjected to different environments only for a fraction of the life-cycle 39 suggests that phenotypic divergent selection during allopatry was intense 11,54 . This finding adds to the knowledge and understanding of adaptive variation in VN among fish populations, a phenomenon that previously was documented mainly for the Atlantic Silverside (Menida menida), for example by Hice, et al. 21 .
To identify the underlying selective driver(s) for the documented divergence would have required a different set of approaches than the ones used in the present study, but we offer two plausible explanations related to the association between VN and swimming performance. Performance studies indicate that VN is associated with locomotion and manoeuvrability 4,31,33,55 , for instance, by influencing body flexibility and the ability to curve the body 56,57 . Moreover, stellar work on sticklebacks proposes that VN mediated differences in swimming performance can translate into differential predation, particularly in juvenile life-stages, such that an optimal vertebral phenotype is favoured 10,34,58 . It is therefore reasonable to hypothesize that dissimilar predation pressures in population-specific breeding habitats may have resulted in divergent selection on VN in pike juveniles. This may ultimately have contributed to subpopulation divergence, in a similar manner as previously suggested for adaptive divergence in growth trajectories of pike 39 . Another explanation is that the subpopulation divergence in VN may have arisen due to the endeavours of the migratory life-history of pike. Given that VN influences swimming performance (as discussed above) it is plausible that VN can also influence migratory success and the ability to cope with physical obstacles faced during migration 33 . The variation in VN among subpopulations may therefore in part reflect differences in physical obstacles, for example rapids, which must be passed to reach breeding sites 55 .
The adaptive significance of individual variation in VN has been emphasized by studies reporting on associations between vertebral phenotypes and fitness correlates 2,4,34,58,59 . Much focus has been on whether VN is positively associated with body size 11,20,59 . We found a curvilinear relationship between VN and body size in juveniles indicating that an intermediate trait value was optimal, similarly to what was shown by Arnold 4 in the garter snake. This pattern might reflect that locomotor performance peaks at intermediate number of vertebrae and that this translates into higher foraging success 10 . Our data further revealed that the relationship between VN and body size was expressed more strongly in the natural environment than in the laboratory, as previously shown in snakes 36 . A reasonable explanation for such context dependence is that any differences between vertebral phenotypes in performance (mobility) would likely be of higher value in a natural more challenging environment, where individuals need to capture free-moving prey, than in the laboratory. An additional explanation is that the inferior performance of individuals with extreme VN might partly reflect correlated responses to environmental stress (see below), and that any malfunctions associated with perturbed developmental processes manifested more strongly in the wild.
Related to the context dependent effects of VN on body size discussed above is the suggestion that the influence of phenotypic trait values on fitness can vary depending on sex, size class or age 3,4,7,34 . Still, few studies have evaluated whether effects of variation in VN on body size persist through time. By reconstructing growth trajectories using cleithra, we demonstrated that the curvilinear relationship between VN and body size documented in juveniles was present at least until individuals reached maturity. It has been suggested, but not empirically tested, that variation in VN also can influence reproductive allocation strategies independently of any effects of VN on body size 9,11,24,34,60 . We tested this hypothesis and found that female reproductive effort adjusted for body size (GSI) reached a peak in individuals with an intermediate number of vertebrae. It remains to be investigated whether this applies also to male reproductive investment.
The signatures of stabilizing selection on body size, growth rate and reproductive investment discussed above might reflect a combination of the aforementioned direct effect of VN on mobility and performance 4,10,31-33 and indirect effects stemming from influences of environmental stress on gene expression and developmental processes 61 . Available evidence indicates that relatively high and relatively low VN are associated with exposure to extreme temperature conditions during early embryonic development 4,18,24 . If extreme temperatures modify gene regulation, activate otherwise unexpressed genes, and impair canalization this can induce developmental perturbations and lead to fitness reducing phenotypic abnormalities 61 that may be related to but not caused by VN as such.
Given the associations with mobility and body size (growth rate), VN is likely subjected to viability selection 4,10 . We compared juvenile and adult life-stages and found that VN decreased across life-stages in all three subpopulations, indicative of directional viability selection. The duration of sampling ranged between 1 and 4 years depending on population and life-stage, and our samples of adults included individuals of different age (∼ 3-10 years). It is therefore unlikely that extreme environmental conditions during certain year(s) systematically affected VN distributions, or that the comparisons of different cohorts confounded environmental influences on vertebral development with selection. Our result also is congruent with previous studies on three-spine sticklebacks 10,58,60 .
Taken together, our longitudinal and cross-sectional analyses of associations between VN and different fitness components suggested that the variation in VN within populations is influenced by a combination of stabilizing and directional selection. That the mode of selection differs among fitness components is common 34,62-64 , and may reflect trade-offs among different components of fitness, such that trait values that enhance survival may impair reproductive effort or growth 3,34,50,65,66 . Moreover, previous studies (e.g. 34,67 ) show that the optimal VN phenotype shifts with changes in body size and/or age. This would result in stabilising selection within life-stages (i.e. adult or juvenile) but directional selection between life-stages, in agreement with our findings.
In conclusion, results showed that differences in VN in pike reflected a combination of genetic and environmental (temperature) effects, and influenced aspects of individual performance (body size, growth rate, reproductive effort and survival) that likely contribute to variation in lifetime reproductive success. The overall patterns could be accounted for by a combination of divergent selection contributing to differentiation among subpopulations, and stabilizing selection reducing variation within subpopulations. Future studies should investigate the contribution of any correlated and indirect effects associated with perturbed development to patterns of variation in VN and fitness components in fish and other organisms. It is noteworthy that signatures of adaptive divergence (which were of relatively large magnitude compared to previous studies conducted within similar spatial scales) were evident despite that subpopulations were only separated by short distances relative to dispersal capacity for a brief period during reproduction and early development 39 . Collectively, our findings expand the present understanding of adaptive variation at fine spatiotemporal scales. In an era of accelerating anthropogenic influence on ecosystems such knowledge is key to successful conservation and management of biodiversity and ecosystem functioning.

Methods
Study system and sampling procedures. Pike is a suggested keystone top-predator fish that occupies a wide variety of aquatic habitats in the northern hemisphere 37 . Typically pike forage using an ambush strategy 43 and it is prone to attack prey that can be larger than half their own body length 37 . Pike has indeterminate growth and breeds yearly after reaching maturity at age 1-4 37 . In some systems, for example the Baltic Sea, pike undertake breeding migrations in excess of 10 km to reach suitable breeding habitats and reproduction occurs in shallow waters during spring (March-May) at water temperatures between 6 and 14 °C 37,41,68 .
For this study, we used six subpopulations of pike breeding in streams entering the Baltic Sea (salinity is around 0.7% in the coastal area) within a confined geographical area with a radius of 25 km (Fig. 1a; Table 1). These streams are known as reproductive and nursery habitats for genetically differentiated subpopulations of anadromous pike 40  Swedish mainland. Kårehamnsbäcken (I) was equally small as streams L, O, T and D but flowed through the agricultural landscape of the island of Öland (Fig. 1a). Marking studies to monitor timing and direction of migration have indicated that the majority of adult pike originating from these streams have a migratory life-history where they only entered the streams during breeding periods and subsequently emigrated back to the Baltic Sea within a few weeks 40,41,69 .
We sampled juveniles and adults in breeding streams to study patterns, causes and consequences of variation in the number of vertebrae, both within and among subpopulations. Adults (N = 328) were captured using fyke-nets placed in the streams during breeding migration (March-May). Individuals were measured for total body length (to the closest cm), from tip of snout to end of caudal fin, using a measuring board. Juveniles (less than six weeks old) were captured (N = 126, Table 1) during emigration from the streams to the Baltic Sea using either hand-trawl or larval traps. Juveniles were measured for total body length using digital callipers (to the closest mm).

Methods for counting the number of vertebrae.
We counted the number of vertebrae in 721 (328 adults and 393 juveniles) individuals originating from six subpopulations and representing both wild-caught and laboratory reared individuals ( Fig. 1a; Table 1). In adults, we counted the vertebrae by analysing radiographs (generated by the x-ray system Arcoma PanoRad, Mediel accessed at the County Hospital of Kalmar, Sweden) in the software program ImageJ and/or by dissection after all soft tissues were removed by boiling (Fig. 1b). In juveniles, vertebrae were visualized by staining the vertebral column using Alizarin red S (Acros Organics, New jersey, USA) according to a method described by Connolly and Yelick 70 followed by dissection. Specimens were then photographed in a dissection microscope (Olympus Digital camera UC 30, Olympus Microscope SZX12) and vertebrae were counted by image analysis using ImageJ software (Fig. 1c). To determine the repeatability of VN counts, 69 individuals (N adults = 39, N juveniles = 30) were counted twice and blind with respect to the first count. Measurement repeatability (intraclass correlation coefficient) 71 was very high for both life stages (adults: 99.1%, F 38,39 = 114.16, P < 0.001; juveniles: 96.7%, F 29,30 = 31.71, P < 0.001), with less than 1% and 4% of the variance in VN being due to measurement error in adults and juveniles, respectively.  Fig. 3). Data were analysed by Fisher's Exact test implemented in SPSS 22, and performed separately for adults and juveniles. This was done to avoid any confounding effects of differences in VN between life-stages (see Results). To avoid cells with small expected values, data of vertebral number were grouped as ≤ 59, 60, 61, or ≥ 62 34 .

Disentangling sources of variation among subpopulations.
To address whether genetic modifications accounted for variation in VN among subpopulations, we compared the frequency distribution of VN among juveniles reared in a common garden experiment. For detailed methods of this section see 39 . In short, we dry-stripped gametes from spermiating males (N = 35) and ovulating females (N = 33) that were captured (within 11 days to reduce any effects of reproductive timing) during the breeding migration in three streams (subpopulations L, O & A; Table 1). Gametes were put on ice and transported to the laboratory where gametes were dry-mixed and artificially fertilized by adding water originating from the experimental facility. Temperature (14 °C) and light regime (14L/10D) were held constant throughout the experiment. To estimate the contribution of additive genetic variation to differences in VN, we used a maternal half-sib breeding design where each female was crossed with two males. Each unique male-female crossing contributed with five aquaria replicates containing five larvae each. Fertilized eggs hatched after 5-6 days of incubation, and offspring were then reared for 35 days and fed brine shrimp (Artemia salina) and cyprinid larvae (Cyprinius sp.) ad lib before the experiment was terminated and juveniles were measured for total body length (TL, mm). For this study, we counted VN in a random subsample of individuals produced in the common garden experiment (one aquaria replicate from each parental crossing, N = 203 juveniles).
To analyse whether frequency distributions of VN differed among subpopulations when reared in a common garden, we used Fisher's Exact tests due to low expected frequencies and data of VN were grouped as presented above. We also compared the frequency distribution of VN in the common garden experiment to that in wild-caught juveniles to infer the plastic component by performing Fisher's Exact test separately for each subpopulation. Analyses were performed in SPSS 22. Testing for developmental plasticity of vertebral number in response to temperature during embryonic development. To examine whether the variation in VN among individuals and subpopulations could be accounted for by differences in temperatures experienced during early embryonic development we performed a split-brood laboratory experiment and incubated eggs under four different temperature regimes (5, 10, 15 and 20 °C). This experiment was designed to also test for gene by environment interaction effects (G × E) on VN, indicative of genetic variation in reaction norms required for evolution of phenotypic plasticity. The experiment was conducted by first dry-stripping gametes from five females and five males originating from subpopulation O (Table 1). Gametes were crossed to obtain five full-sib family batches. Family egg batches were artificially fertilized (conducted according to the procedure described for the common garden experiment), split into four parts consisting of ca. 20 eggs each, and incubated in 2L aquaria in four temperature treatments with each family being represented in each temperature. After hatching, larvae fed on the yolk-sac until depletion and were thereafter provided brine shrimp (Artemia salina) ad lib. Individuals were reared until they reached a length of more than 25 mm before termination to facilitate counting the number of vertebrae after staining with Alizarin red S and dissection. The experimental temperature gradient (5-20 °C) was chosen to comprise the general range of breeding temperatures for pike 37,68 . Data on VN were analysed using procedure GLM in SPSS 22 by including temperature treatment, family and the treatment by family interaction as explanatory fixed factors.

Evaluation of whether subpopulation divergence in vertebral number is adaptive by Q ST -F ST comparison.
To assess whether divergence among subpopulations in VN could be more readily accounted for by differences in selection pressures than by random genetic drift, we compared phenotypic differentiation from the common garden experiment (Q ST ) with genetic differentiation at neutral microsatellite loci (F ST ). For this comparison we used data for populations L, O and A (Table 1). Q ST was estimated according to the method described by Spitze 46 (i.e. Q ST = σ 2 b /(2σ 2 w + σ 2 b ), where σ 2 b and σ 2 w are additive genetic variances between and within subpopulations 47 ). Restricted maximum likelihood (REML) estimates of σ 2 b and σ 2 w were derived from a Linear Mixed model applied to the common garden data using procedure MIXED in SPSS 22. In this model, VN was set as the dependent variable, subpopulation as fixed factor, and female and male as random factors nested within subpopulation. To avoid confounding effects of non-genetic variance, only variance components based on males were used to estimate Q ST . The 95% CIs for Q ST were estimated using the direct simulation data method 72 . The estimated neutral genetic differentiation (F ST ) among subpopulations used in the common garden experiment was 0.069 (Bootstrapped 95% confidence interval = 0.048-0.085) as presented by Tibblin, et al. 39 .
Testing for associations of variation in vertebral number with fitness components. We examined whether VN was associated with survival by a cross-sectional comparison of the frequency distribution of VN between wild-caught adults and juveniles representing three subpopulations (L, O & D, Table 1). Data were analysed using a Log-linear model as implemented in procedure GENLIN in SPSS 22 to be able to control statistically for any variation in vertebral number among subpopulations. Vertebral number (grouped as ≤ 59, 60, 61, or ≥ 62) was treated as dependent variable with a Poisson distribution, and life-stage (i.e. adults or juveniles), subpopulation and the life-stage by subpopulation interaction as explanatory factors. The likelihood ratio was set as method for the Chi Square statistics. The interaction was non-significant (P = 0.30) and subsequently deleted which resulted in a final model without interaction.
We evaluated whether and how variation among individuals in VN was associated with juvenile growth rate by analysing the relationship between body size (mm length) and VN in juveniles that had been reared in captivity (from the common garden experiment described above), and in juveniles that were captured in the wild upon leaving their breeding stream. To allow for comparisons between VN and different fitness components (see below), we used data only for juveniles originating from subpopulation O, which is the population with the largest sample size. Data were analysed using multiple regression as implemented with procedure REG in SAS. Body size was set as response variable, and VN and squared VN was included in the model to test for linear and curvilinear (quadratic regression) relationships with body size 48 . We first analysed captive-reared and wild-caught juveniles separately. Next, we evaluated the hypothesis that effects of VN on performance and growth might be context dependent and expressed more strongly in more challenging environments where manoeuvrability is more important. To this end, we analysed data for captive-reared and wild-caught juveniles simultaneously and tested for effects on body size of the interaction between origin and VN using procedure GLM in SAS. To further evaluate the evidence for a curvilinear effect we performed a randomization test procedure for quadratic regression. The observed values of body size were randomly reallocated to the VN values, a 2-tailed test for a non-zero quadratic effect was conducted, and the outcome of 2000 randomizations was used to compute an overall significance level using the SAS macro RTEST 73 .
Scientific RepoRts | 6:26372 | DOI: 10.1038/srep26372 To evaluate whether and how VN influenced individual growth rates throughout the first three years of life (sub-adult life-stage) we reconstructed growth trajectories by measuring annual growth rings (zonation patterns) in the cleithrum for a subsample of adults originating from subpopulation O (N = 108). We restricted the analyses to the first three growth rings to minimize confounding effects of any differences in reproductive strategies. In pike, length of cleithrum is linearly related to total body length and constitutes a reliable proxy of growth rate (for details see 39 ). To test for associations between individual growth trajectories and VN we first calculated annual growth rate for each individual, computed average growth rate for all individuals that had the same VN, and then performed a multiple regression analysis based on the mean values to test for linear and curvilinear (quadratic) effects.
To evaluate whether and how VN influenced reproductive investment, we first calculated GSI (Gonad Somatic Index) for a subsample of individuals originating from subpopulation O (N = 48, Table 1). We sampled individuals when they were entering the spawning areas and only used ripe but not ovulating females. We are therefore confident that we sampled individuals at the peak of their GSI 37 . To calculate GSI we divided gonad wet mass with somatic body mass (measured with 10 g accuracy, Digital Scale E15302, Fox International Group Limited) 49 . Next, we estimated least squares means of GSI for all individuals that had the same VN. These were obtained from an ANCOVA model in which GSI was treated as the dependent variable, female body length was treated as a covariate, and vertebral number and year were treated as class variables (the model was significant: F 9,37 = 2.28, P = 0.038). We then performed a multiple regression analysis based on the least square means to test for linear and curvilinear (quadratic) effects 4,31 .
Ethics statement. Ethical approval for the study was granted by the Ethical Committee on Animal Research in Linköping, Sweden (approval 9-06 & 39-10). The methods were carried out in accordance with the approved guidelines.