Striking Phenotypic Variation yet Low Genetic Differentiation in Sympatric Lake Trout (Salvelinus namaycush).

The study of population differentiation in the context of ecological speciation is commonly assessed using populations with obvious discreteness. Fewer studies have examined diversifying populations with occasional adaptive variation and minor reproductive isolation, so factors impeding or facilitating the progress of early stage differentiation are less understood. We detected non-random genetic structuring in lake trout (Salvelinus namaycush) inhabiting a large, pristine, postglacial lake (Mistassini Lake, Canada), with up to five discernible genetic clusters having distinctions in body shape, size, colouration and head shape. However, genetic differentiation was low (FST = 0.017) and genetic clustering was largely incongruent between several population- and individual-based clustering approaches. Genotype- and phenotype-environment associations with spatial habitat, depth and fish community structure (competitors and prey) were either inconsistent or weak. Striking morphological variation was often more continuous within than among defined genetic clusters. Low genetic differentiation was a consequence of relatively high contemporary gene flow despite large effective population sizes, not migration-drift disequilibrium. Our results suggest a highly plastic propensity for occupying multiple habitat niches in lake trout and a low cost of morphological plasticity, which may constrain the speed and extent of adaptive divergence. We discuss how factors relating to niche conservatism in this species may also influence how plasticity affects adaptive divergence, even where ample ecological opportunity apparently exists.

Introduction

Ecological speciation occurs when divergent natural selection in different environments results in reproductive isolation and ultimately new species formation [1]. Its speed and extent across taxa likely depends on multiple factors meriting further investigation [1,2]. Hendry [2] suggested that ecological speciation is best considered as a ‘speciation continuum’, wherein species range from adaptive variation between randomly mating populations all the way to complete and irreversible reproductive isolation between two distinct species. A better understanding of this continuum in evolutionarily young populations of harvested species could improve biodiversity conservation, because evolving phenotypic and genetic differentiation can significantly influence population growth, productivity and species persistence [3,4].

Many north temperate freshwater fishes occupying postglacial lakes, including harvested species, contain phenotypically- and genetically-distinct morphs that have primarily diverged in trophic niche use [2,3,5-7]. They hence provide excellent models for both the study of differentiation along the ecological speciation continuum and its consequences for biodiversity conservation. In several species, such as Arctic charr (Salvelinus alpinus) and threespine stickleback (Gasterosteus spp.), morphs are at the intermediary stages of ecological speciation. These morphs exhibit varying degrees of genetic differentiation, moderate to relatively high gene flow, adaptive variation to limnetic and benthic feeding niches, and potentially reversible reproductive isolation [2,8-10]. In more extreme cases such as lake whitefish (Coregonus clupeaformis), divergent morphs exhibit complete or near complete reproductive isolation, and genetic incompatibilities might arise as a result of divergent selection, or for other reasons due to divergent selection following gene flow reduction [2,11].

Less empirical attention has been directed to diversifying postglacial fishes that exhibit occasional adaptive variation with minor reproductive isolation, and hence to the potential factors influencing the progress or lack of progress at earlier stages of differentiation (but see [10,12]). One species that might typify this situation is lake trout (Salvelinus namaycush, Walbaum, 1792), a predator species which exhibits sympatric morphs associated with occupying different habitat niches within several, large North American postglacial lakes [13,14]. The best known morphs are the ‘leans’, ‘siscowets’ and ‘humpers’ of the Laurentian Great Lakes, which occupy mid-depth, deep-water, and drop-off habitats, respectively. The morphs have trophic specializations for consuming different prey [14-16], but appear to exhibit weak genetic differentiation [17]. Morphs in a few other lakes have similarities and differences to Great Lakes morphs. In Great Bear Lake, morphological variation between piscivorous and insectivorous morphs is unrelated to depth [13,18] and up to four weakly, genetically differentiated morphs inhabit shallow waters [19,20]. Conversely, in Atlin Lake, two morphs are not genetically-distinct but have clear depth and habitat preferences [21]. Finally, in Great Slave Lake, three identified morphs display trophic specialization to benthic and pelagic habitats but their genetic differentiation has yet to be studied [14,22]. Overall, previous research suggests that lake trout may be at the earlier-stage end of the ecological speciation continuum in multiple lakes.

Mistassini Lake (50.958034°, -73.622340°; Fig 1A) is a large (2,335 km2), deep (180 m), postglacial lake in northern Canada [23]. It supports numerous habitat niches that has favoured sympatric differentiation in several harvested fish species, including brook trout (Salvelinus fontinalis) and walleye (Sander vitreus) [24,25]. For lake trout, Wilson & Hebert [26] suggested that multiple mitochondrial DNA lineages colonized Mistassini Lake after the last deglaciation (7,000–8,000 BP). Two morphs differing in growth, age- and size-at-maturity and depth use were also recently described [27,28]. While these studies offered a glimpse into what such a large and pristine lake might harbour, the small sample sizes, limited spatial distribution and limited genetic assays left an incomplete picture of sympatric differentiation.

Fig 1

Geographic location of Mistassini Lake, Québec, Canada and colour frequency, depth distribution and basin preferences of each one of the genetically-demarcated clusters.

(A) Shown on the map are 24 hr gillnet sets (star) and lake sector (W1–6 and E7–12) divided by dashed lines. (B) Frequency of colours observed in each genetically-demarcated lake trout clusters in Mistassini Lake. (C) Depth distribution of lake trout in each one of the genetically-demarcated clusters. Beanplot shows accurate densities and the red line indicates the median. (D) Basin preferences of lake trout in each one of the genetically-demarcated clusters. Results of a generalized linear model (GLM) that took into account the fishing effort and demonstrates that clusters 1 and 5 were caught in higher numbers in the eastern basin and clusters 3 and 4 were caught in higher numbers in the western.

We firstly tested and provided evidence that Mistassini Lake harboured multiple genetically- and morphologically-distinguishable clusters of lake trout associated with spatial habitat, depth and fish community structure, through the use of morphological analyses, catch locations, bycatch data and microsatellite DNA surveys. Yet genetic differentiation was low and genetic clustering was largely incongruent among several population- and individual-based clustering approaches. Moreover, genotype- and phenotype-environment associations were either inconsistent or weak. We thus also evaluated genetic differentiation in relation to effective population sizes and the extent (and possible asymmetry) of contemporary gene flow, and considered additional ecological factors within Mistassini and other large postglacial lakes that might inhibit genetic differentiation in lake trout and other postglacial fishes. Our work contributes to a better understanding of how intraspecific phenotypic plasticity might affect the trajectory and progress towards ecological speciation.

Materials and Methods

Fish sampling

Between June 12 and July 29, 2013, 636 lake trout were sampled throughout Mistassini Lake using monofilament nylon gillnets and angling, in order to obtain genetic (n = 636) and morphological data (n = 281). Gillnet design was based on previous lake trout research in large lakes (183 m long, 1.8 m tall; six 30.5 m gangs, one for each mesh size, ranging 51–114 mm stretch mesh; [27], and permitted unbiased targeting of all known sizes of Mistassini lake trout [29]. Twenty-nine gillnets were set throughout the lake (Fig 1A) in areas known by local fishers to have lake trout, from depths between 3–178 m, were soaked overnight and lifted after a 24-hour period. Gillnets were lifted slowly, bycatch were recorded and returned to the lake immediately, and trout were placed in fresh water baths with aerators. All living trout were anaesthetized with buffered tricaine methanesulfonate (MS-222) and subsequently processed. From each trout, we collected and preserved a small piece of adipose fin tissue in 95% ethanol until DNA extraction, took a standardized photograph (see below), measured total and fork length (TL and FL, respectively, ± 1 mm) and wet mass (± 50 g), and recorded net or lure depth and location of capture (GPS). Tissue samples were also donated by local anglers. To account for relative spatial location when GPS data were unavailable for angled samples, the location was reported as a sector, pre-determined on a map provided in a sampling kit. Each basin was separated into six sectors (Fig 1A). Research undertaken in this study complies with the requirements of the Canadian Council on Animal Care (CCAC). This protocol was approved by Concordia University Office of Research–Research Ethics and Compliance Unit (permit number: 30003281) and the Ministère des Ressources naturelles direction de l’expertise Énergie-Faune-Forêts-Mines-Territoire du Nord-du-Québec (permit number: 2013-04-03-107-10-S-P).

Genotyping

We extracted and analyzed DNA from 636 trout by amplifying 19 microsatellite loci using multiplexed polymerase chain reactions (PCRs; see S1 File for loci and PCR details). Amplified products were electrophoretically migrated and allele sizes scored using a 3500x Genetic Analyzer, associated size standards and software (Applied Biosystems Inc.). To estimate repeatability, 16 samples were independently genotyped and scored three times.

Genetic diversity

We tested for the presence of null alleles, large allelic size dropout, or scoring errors associated with allelic stuttering issues in our dataset using MICROCHECKER (v.2.2.3) [30]. Descriptive genetic statistics, including number of alleles (NA), expected (HE) and observed (HO) heterozygosity for each locus were computed using FSTAT (v.2.9.3.2.) [31]. To ensure that loci fulfilled assumptions of selective neutrality in Hardy-Weinberg equilibrium (HWE) tests and in clustering analyses below, tests of selection were performed using ARLEQUIN (v.3.5.1.2) [32]. We then tested for deviations from HWE and linkage disequilibrium (LD) using GENEPOP (v.4.2) [33] (i) under the null hypothesis of one randomly mating population within Mistassini Lake, but also (ii) under the alternative hypothesis that genetic structuring was present, defined as multiple, individual based clusters inferred from the Bayesian model-based clustering software STRUCTURE (v.2.3.4) [34]. If the null hypothesis was rejected, we expected many more HWE deviations to exist when assessing the lake as one cluster vs. multiple ones.

Genetic clustering/differentiation

Non-spatial clustering

We employed both non-spatial and spatially informed (see below) approaches to characterize genetic clustering and structure. Non-spatial analyses were performed using three different methodologies (see S2 File); however, we only report the methods and results from STRUCTURE in the main paper, as the most commonly used clustering software with more easily biologically-interpretable results. Namely, this software consistently defines clusters (K) [35] and accounts for admixture [34], which is important for situations of low FST and multiple glacial lineages. We ran STRUCTURE with no a priori information under models of K = 1–20, a burn-in period of 500,000 followed by 2,000,000 iterations, and replicated 20 times each per K. We estimated the most likely K value based on the ad hoc ΔK statistic of Evanno et al. [35] as implemented in STRUCTURE HARVESTER [36] and the highest mean value of the log-likelihood of the data. The 20 iterations per K value were then combined into a single output, providing an overall inferred ancestry coefficient (q) for each individual using the program CLUMPP (v.1.1.2) [37]. Because STRUCTURE can over estimate K when closely related individuals are present in a sample [38] we estimated the number of full siblings using COLONY (v.2.0.6.1) [39] as a supplementary analysis, assuming male and female polygamy and default parameter settings over three runs. We then randomly removed one half of each full sibling pair and reran STRUCTURE with the same parameter settings as described above. COLONY results revealed a low number of full sibling pairs (8, i.e., 16 of 636 individuals) from separate families and their exclusion did not affect the mean LnP[D] output from STRUCTURE (data not shown). Therefore, all subsequent analyses were based on the full dataset. The extent of genetic differentiation (FST) [40] between defined clusters was quantified using GENEPOP. We also calculated pairwise FST for both body and head morphological clusters (see below) to determine if genetic differentiation was higher among morphological vs. genetic clustering.

Genetic differentiation may be influenced by stepwise mutations that are accounted for in measuring RST rather than FST, which in young postglacial lakes can reflect population divergence predating lake age [41]. We assessed and confirmed that RST > FST in Mistassini Lake (S3 File). Nonetheless, to make our study comparable to previous lake trout research, subsequent analyses were based on FST only, because FST (i) typically has lower standard errors than RST [42], (ii) is more precise when high gene flow is detected [43], and (iii) performs better than RST for most typical sample sizes [44].

Spatial differentiation

To complement non-spatial clustering analyses and provide insight into individual-level genetic relationships, we tested for associations between geography (longitude, latitude and depth) and the distribution of individual genetic variation across the lake. The advantage of these spatial analyses is that they do not rely on HWE or linkage equilibrium and thus are useful for detecting gradients and complex spatial-genetic structures [45]. We were able to obtain GPS (± 5 m) and depth (± 1.8 m) on 554 of 636 trout for these analyses. Similarly to above, we performed spatially-informed analyses using three different methods (see S2 File); however, we only report the methods and results from a spatial principal components analysis (sPCA) performed using the ‘adegenet’ package (v.1.3–1) [46] in R [47]. An sPCA compares the allelic frequency of individuals to allelic frequencies of their neighbours using Moran’s index (I) [45,48]. Although our samples were collected relatively evenly across the lake, some were more heavily sampled and so we used two different connection networks to define which individuals were neighbours in the sPCA algorithm: (i) Delaunay triangulation [49], which is suited for uniform sampling and commonly used in a number of different spatial programs [50]; and (ii) inverse distances, which is recommended when sampling is aggregated [45] and assumes all individuals are neighbours. We tested for significant global (i.e., neighbouring individuals are more similar than expected) and local (i.e., neighbouring individuals are more dissimilar than expected) structures using 999 permutations.

Contemporary effective population sizes and gene flow

We estimated the contemporary effective population size (Ne) for each defined genetic cluster using the linkage disequilibrium method implemented in LDNe (v. 1.31) [51]. Lake trout were not aged, so our cluster samples consisted of overlapping generations which downwardly biases Ne estimates up to 25–30% [52]. We estimated contemporary gene flow (recent migration rates, m, within the last few generations) between cluster pairs using all individuals with BayesAss (v.3.0) [53]. Five separate iterations were performed using 107 Markov chain Monte Carlo (MCMC) iterations with a burn-in period of 106 and sampled at a frequency of 100. Prior to each run, mixing parameters (allele frequencies, inbreeding and migration rates) were optimized until acceptance rates were between 20–35%, as recommended in the manual. To assess whether MCMC had converged successfully for each run and overall, we used the program TRACER (v.1.6.0) [54]. Whether mutation-drift equilibrium conditions had been reached among defined genetic clusters was assessed using Whitlock’s [55] equation for the time, in generations, required for FST to reach halfway to a new equilibrium:

, where m = mean gene flow into each cluster.

Morphological analysis

Body and Head Morphology

To characterize body morphological variation, we took a full-bodied standardized photograph of each trout immediately after capture using a digital camera (Nikon D3100) with a UV filter mounted on a tripod. Many tissue samples were donated from local fishers, so we were only able to take photos of a total of 281 of 636 trout. We placed each fish on its right side on a flat piece of plywood with the dorsal, caudal and anal fins in open positions; if fish were bent or distorted, photographs were discarded. Eighteen landmarks (Fig S4.1A in S4 File) were then digitized on each trout using tpsDig2 (v.2.17) [56]. A separate morphological analysis was performed on the head with three landmarks and 22 semi-landmarks used to measure curvature, following Zimmerman et al. [57] (Fig S4.1B in S4 File). Ten equally-spaced regions were produced by a reference grid and superimposed on each photograph between the snout and opercle using MakeFan (v.8) [58]. Semi-landmarks were then slid along both upper and lower curves using tpsDig2, aligning semi-landmarks perpendicularly to the curve by reducing the bending energy among individual points [59]. Landmark and semi-landmark positions utilized are common in studies of lake trout and related species, and were selected to measure traits commonly associated with adaptation in swimming and foraging performance [8,24,57,60].

Geometric morphometric analyses were conducted separately on body and head morphology datasets using tpsRelw (v.1.36) [61]. TpsRelw takes into account spatial variation among assigned landmarks relative to all others using a thin-spline analysis [62]. It specifically produces partial warps, which are geometric constructs derived from the amount and direction of bending required to change the consensus shape [59]; principal components of these partial warps scores are termed relative warps (RWs), which quantify body shape and head variation [63]. The first seven RWs for body shape (78% of variation) and first three RWs for head shape (70% of variation) were used for subsequent statistical analyses. We identified morphological clusters with MCLUST (v.4.4) [64] as implemented in R. Tested models assumed that the data represented one to up to 10 clusters. All models were analyzed using two (EII, VII) of the nine multivariate mixture options; the model with the highest Bayesian information criteria (BIC) was selected as best describing morphological clusters [65]. To measure the reliability of classification of individual fish to MCLUST cluster, we performed body and head discriminant function analyses (DFA) with jackknifing in PAST (v.3.08) [66].

Body size and colouration

Body size (mass) and the relationship between mass and total length (TL) was compared among genetic clusters using linear models. In the latter model, the natural logarithm (ln) of mass was regressed onto the ln of TL, genetic cluster and their interaction. A subsequent pairwise comparison of the inferred slopes (mass/length curves) was conducted using ‘lsmeans’ package (v.2.17) [67]. Mass was recorded for 520 and TL for 472 (of 636) trout.

Colour frequency distributions were compared using a χ2 test in the ‘vcd’ package (v.1.4–1) [68]. The colour of 422 (of 636) trout was classified and recorded at the time of capture based on being black, brown, dark silver, silver or light (i.e., white, light silver).

Associations between genetic, morphological and ecological differentiation

Genetic & morphological or ecological associations

We tested for an association between genetic and morphological clusters using a contingency test. We then assessed whether the spatial distribution of genetic clusters was influenced by geography (lake basin, sector within basin) and/or by water depth, using a redundancy analysis (RDA) within the ‘vegan’ package (v.2.3–1) [69], which incorporated raw q values for each individual. We used a global permutation test of the RDA result to test for a relationship between genetic cluster and ecological variables, as well as a permutation test on canonical axes with 1,000 steps [70]; R2 values were calculated as a measure of ‘goodness-of-fit’. The RDA was complemented with a similar multinomial regression analysis which, in contrast, incorporated all individuals weighted based on their highest q value. Because the RDA and multinomial regression analyses yielded largely congruent results, only results for the former are reported in the paper (see S5 File for multinomial regression details). Differences in the spatial distributions of clusters among basins were also assessed using a generalized linear model (GLM) fitted with a binomial error distribution and logit link-function (proportion data).

Genetic & fish community structure associations

Based on recorded gillnet bycatch data, we tested whether genetic clusters were associated with fish community structure using a contingency test. Bycatch abundance and diversity captured within the same net/panel as each trout provided a reasonable proxy for prey availability and the extent of interspecific competitor species. The relationship between genetic clustering and fish community structure was visualized using a RDA.

Mantel tests: genetic differentiation vs. morphological or ecological differentiation

We further tested for associations between the extent of genetic differentiation among defined clusters and either morphological differentiation, depth, fish community structure (abundance), or colour dissimilarity, by using Mantel tests implemented in the ‘vegan’ R package. Mantel tests always compared two pairwise distance matrices: cluster pairwise FST values and either the absolute difference in (i) mean scores for the first seven RWs for body morphology and first three RWs for head morphology, (ii) median depth, (iii) prey abundance, or in (iv) colour frequencies.

Statistical notes

Unless otherwise stated, we corrected multiple comparison P-values to adjust for any potential type I errors following the false discovery rate procedure [71]; significance was defined at the nominal alpha value of 0.05.

Results

A high repeatability of genotyping (100%) was confirmed for the 16 samples independently run and scored three times. The 19 loci averaged 29.5 alleles per locus (range of 13–72) with an average HO across loci of 0.813 (range of 0.526–0.939; Table S6.1 in S6 File). No locus showed evidence of being under selection, nor was there any evidence of scoring errors as reported from MICROCHECKER. Under the null hypothesis that all 636 trout were from one randomly-mating population, heterozygote deficiencies were detected at 15 of 19 loci (P < 0.050; Table S6.2 in S6 File). These deficiencies were greatly reduced at the level of clusters defined by STRUCTURE, providing a first indication of genetic structuring (28 of 95 cluster-locus comparisons at uncorrected P-values (0.000–0.049); the remaining 17 deviations from HWE (P = 0.000–0.044) were spread across all five clusters and across nine loci. We also detected significant LD in only seven of 171 locus-by-cluster tests (20 of 171 tests at uncorrected P-values), and found little evidence for LD when accounting for clusters (12 of 855 tests; Table S6.3 in S6 File). Significant LD tests were spread across all clusters and among 23 unique loci-pairs.

Results from STRUCTURE and STRUCTURE HARVESTER suggested five genetic clusters in the lake (ΔK = 57.03 mean LnP[D] = -58060.42; Table S6.4 in S6 File), with a global FST of 0.017 and significantly greater than zero. Pairwise estimates of FST between clusters ranged from 0.012 to 0.036 (Table S6.5 in S6 File); cluster 2 vs. 3 and 4 vs. 5 exhibited the lowest and highest levels of differentiation, respectively. Pairwise estimates of FST between morphological clusters (see below) were lower and ranged from 0 to 0.002 (body; Table S6.5 in S6 File) or 0.003 (head)).

The global permutation test on the eigenvalues derived from the sPCA using a Delaunay connection network revealed significant global structure (nper = 999, P = 0.028) but no significant local structure (nper = 999; P = 0.621). Screen- and bar-plots of the eigenvalues (Fig 2) suggested that only the first global axis should be retained. This axis revealed genetic complexes in which individuals were more closely related than expected, with particular nuances existing in the western basin near the mouth of the Rupert River and in the north near the Big Pass (Fig 2). Using the inverse distance connection method did not significantly change the results (not shown).

Fig 2

Spatial principal component analysis of lake trout showing the first global structure across Mistassini Lake.

Plotted are the lagged scores in which colours (blue to red) represent the score of an individual genotype; each is positioned by its spatial coordinates. Inset is the barplot for all eigenvalues (A), and the screen plot (B) which illustrates the spatial and variance components of those eigenvalues.

Point estimates of Ne in general had fairly tight confidence intervals and ranged from 371 (cluster 5) to 4,199 (cluster 2; Table S7.1 in S7 File). Estimates of m were quite high, ranging from 0.002 to 0.091, with a global average of 0.017 (Table S7.2 in S7 File). Asymmetries in m existed, namely higher m was detected from cluster 1 to clusters 2, 4 and 5 than vise-versa (Fig S7.1 in S7 File). Generation times (t) for FST to reach a new equilibrium between cluster-pairs ranged from 22 (clusters 2–5) to 177 generations (clusters 1–3), suggesting migration-drift equilibrium has been reached.

Morphological Analysis

Body and Head morphology

The seven body RWs explained 78% and the first three head RWs explained 70% of the total variation (34, 28 RWs, respectively). Major shape differences are visualized in Fig 3 and included the slope of the snout and lower jaw, dorsal and belly curvature (body RW1; 36%), slope of top of cranium and dorsal side, length and depth of caudal peduncle (body RW2; 17%), body depth, eye size (body RW3; 12%), length of head and snout, eye position, length of upper jaw (head RW1; 36%), head depth and bluntness of snout (head RW2; 23%) and slope of lower and upper jaw (head RW3; 12%).

Fig 3

Association between morphological clusters and genetically-differentiated clusters of lake trout present in Mistassini Lake.

Morphological clusters are identified by symbol shape and ellipses which represent 67% of that cluster’s variation and genetically-differentiated clusters are identified by symbol colour. (A) RW1 vs. 2 and 1 vs. 3 (55%) of body shape variation. Morphological shifts for RW1 (26%) correspond to the slope of the snout, lower jaw and dorsal and belly curvature; RW2 (17%) correspond to slope of top of cranium and dorsal side, length and depth of caudal peduncle; RW3 (12%) correspond to body depth, eye position and size. (B) RW 1 vs. 2 and 1 vs. 3 (71%) of head shape variation. Morphological shifts for RW1 (36%) correspond to length of head and snout, eye position and length of upper jaw; RW2 (23%) correspond to head depth and bluntness of snout; RW3 (12%) correspond to slope of lower and upper jaw. Inset images are visualizations of the shape at the most extreme end of each relative warp.

MCLUST defined three separate body shape clusters (B1-3) with a ΔBIC = 42, and five head shape clusters (H1-5) with a ΔBIC = 6 (Table S4.1 in S4 File; Fig 3), and it recommended carrying forward the first seven body RWs body and three head RWs for further analyses. Consensus shapes of the three body clusters did not differ greatly, however, B1 (n = 144) had a long and steep lower jaw and long and narrow caudal peduncle, B2 (n = 111) had a shorter head with a steep cranium top and a long and narrow caudal peduncle, whereas B3 (n = 26) had the greatest body depth, smaller eyes and a thicker and shorter caudal peduncle (Fig S4.2A in S4 File). Morphological variation was greater among head shape clusters. H1 (n = 19) had the bluntest snout and shortest upper jaw. H2 (n = 45) had the longest head and a less blunt and more streamlined snout. H3 (n = 22) was the most streamlined with the longest upper jaw. H4 (n = 177) was very similar to H2 but with a slightly greater depth and longer upper jaw, and H5 (n = 18) had the shortest and deepest head with eyes located most dorsally (Fig S4.2B in S4 File). The DFA analysis resulted in 85% and 92% correct morphological assignment for body and head clusters, respectively (Table S4.2A,B in S4 File).

Mass varied among genetic clusters (F4,515 = 24.88, P < 0.001). Cluster 4 individuals were larger than all others (all P < 0.001); cluster 3 was larger than clusters 1, 2, and 5 (pairwise t-tests; P < 0.001–0.02; Fig S4.3A in S4 File). Both total length (TL) and the TL*population interaction had a significant effect on mass (F9,437 = 575.2, P < 0.001). The slope for the mass/length relationship of cluster 4 was significantly steeper than cluster 1, 2 and 3 (pairwise t-test; P < 0.003–0.006; Fig S4.3B in S4 File).

Genetic clusters varied in colouration type frequencies ( = 130.91, P < 0.001). Clusters 2, 3, and 5 were predominately (between 73–80%) of black colouration (Fig 1B), whereas clusters 1 and 4 contained between 69 and 72% of lighter (dark silver, silver or light) coloured trout.

Genetic clusters were associated with body and head clusters (body: = 38.1, P < 0.001; head: = 38.1, P < 0.002) but these relationships were best described as weak (body: Cramer’s V = 0.26; head: Cramer’s V = 0.18). Each morphological cluster contained individuals in varying frequencies from all genetic clusters (Fig 3).

The selected RDA model included depth, basin and sector (R2adj = 0.12). A global permutation test revealed relationships between ecological variables and genetic clustering (P = 0.001). A subsequent canonical axes permutation test revealed that only the first three RDA axes explained this relationship (all P = 0.001). RDA1 (50% of the variation) was driven primarily by depth and was especially important for distinguishing cluster 1 as occupying deeper water than all other clusters (Fig S8.1 in S8 File). RDA2 (37%) was primarily driven by individuals located in the eastern basin and sectors, and distinguished primarily cluster 5. RDA3 (12%) was primarily driven by individuals located in the western basin and sectors and was most important for distinguishing clusters 2 and 4. A supplementary, pairwise comparison of the least square mean depth showed that cluster 1 was captured at greater depths compared to all other clusters (Fig 1C). Fishing effort-corrected spatial frequencies also exemplified that clusters 1 and 5 were captured disproportionately more in the eastern than western basin (GLM, P < 0.001); clusters 3 and 4 showed the converse pattern (both P < 0.001), whereas cluster 2 was captured in equal numbers among basins (P = 0.886; Fig 1D).

As depth influenced genetic clustering, we also examined whether it affected morphological clustering. Depth was only a weak predictor of body shape cluster (linear model, F2,202 = 3.5, P = 0.033); only individuals in cluster B1 were captured at a greater mean depth compared to B2 individuals (pairwise t-tests, P = 0.027; Fig S8.2A in S8 File). Conversely, depth was a good predictor of head shape cluster (linear model, F4,200 = 7.4, P < 0.001). Cluster H5 individuals were captured at a greater mean depth than all other morphological clusters (pairwise t-tests, all P < 0.025; Fig S8.2B in S8 File); cluster H4 was captured more deeply than clusters H2 and H3 (P = 0.030 and 0.011 respectively) and cluster H1 was captured more deeply than cluster H3 (P = 0.029).

Genetic clusters were associated with different quantities of bycatch of other fish species (P < 0.001), though the association was weak (Cramer’s V = 0.168). Cluster 1 was rarely captured with other fish species, clusters 2 and 3 were associated with burbot (Lota lota) and Lake whitefish, cluster 5 with walleye, and cluster 4 with all bycatch species (Fig S8.3 in S8 File).

There was little evidence for a relationship between cluster pairwise FST and the extent of morphological, water depth and prey differences (Table S8.1 in S8 File). Relationships only existed between FST and the differences in the observed colour frequencies of brown and black (P = 0.008 and 0.017 respectively). Specifically, pairwise FST was positively associated with both dark colours in which clusters with the greatest genetic differentiation also exhibited the greatest absolute difference in the frequencies of individuals that were black and brown. Additionally, some evident positive trends between increasing FST and increasing morphological differences were detected (body: RW1, RW2, RW5, RW6; head: RW1, RW2; Fig 4; Fig S8.4 in S8 File).

Fig 4

Visualization of Mantel tests: genetic distance vs. morphological and ecological variables.

Shown are comparisons of genetic distance (pairwise FST) vs. (i) the absolute difference in mean relative warp (RW) score for the first two RWs for body (A, B) and head morphology (C, D); (ii) absolute difference in prey abundance (E); and (iii) the absolute difference in median depth (m) (F). Body RW1 and 2 represent 26% and 17% of the total variation; whereas head RW1 and 2 represent 36% and 23% of the total variation. The remaining visualizations of Mantel tests can be seen in Fig S8.4 in S8 File.

Discussion

We detected striking phenotypic variation yet low genetic differentiation in lake trout occupying a large, pristine, postglacial lake. Genetic clustering was (i) influenced primarily by depth and less so by lake geography (basin, sector); (ii) associated with some variation in body size, colouration, and fish community structure, but was (iii) inconsistently associated with body/head shape differences linked to trophic ecology and locomotor mode. Rather, genetic clusters displayed high morphological variation, and all contained individuals belonging to the here-identified three body shape clusters or five head shape clusters.

Associations between phenotypic variation, genetic differentiation and the environment; comparisons with other large, postglacial lakes

One obvious phenotype-environment association detected was camouflage body colouration and depth. In contrast to other, clear postglacial lakes inhabited by lake trout morphs, Mistassini is tea-coloured with a relatively shallow secchi depth (~ 9 m; [72]). Correspondingly, genetic clusters associated with shallower waters (#2, 3, 5; median = 5–10 m), mid-depths (#4; median = 27 m) and deeper areas (#1; median = 30 m) were, respectively, predominately black and brown in colouration, a mix of all colours, and frequently silver or white. Second, clusters with the largest (#4) and smallest (#1) body sizes were most frequently captured, respectively, at mid-depths associated with steep slopes harbouring all prey species and with less abundant prey species.

Functionality of different head shapes also corresponded somewhat with depth. Lake trout with dorsally located eyes, blunt snouts and short heads (H5) were routinely captured in deep water (median = 142 m), whereas individuals with more streamlined heads (H1 –H4) were captured in shallower water (median = 9–27 m); the former head shape is specialized for benthic foraging or more vertical oriented piscivory, the latter is associated with piscivory or insectivorous preferences in shallower waters [9]. Body shape differences between morphological clusters were less obvious and characterized by changes in body depth and caudal peduncle dimensions. Individuals captured in deeper water (B1 and B2; median = 33 and 45 m respectively) had long and narrow caudal peduncles which is needed for sustained, pelagic swimming [73]. Conversely, trout captured in shallower water (B3; median = 31 m) had the greatest body depth and much thicker and shorter caudal peduncles—traits usually linked to enhanced manoeuverability in more complex habitat [74], burst swimming and acceleration [73].

These general body shape-depth/habitat associations have similarities and differences with other lake trout morphs. The shallow-water, streamlined, long and narrow caudal peduncle ‘lean’ morph in other lakes [14] is similar to some of Mistassini’s B1 and B2 (and shallow-water genetic clusters 2, 3 and 5), and likely derived for pelagic swimming. Two deepwater morphs identified in other lakes–small-bodied ‘humpers’ with large, dorsally positioned eyes and typically found on steep slopes, and deeper-bodied ‘siscowets’ that have a sloping snout and thick and short caudal peduncles [14]–have shared similarities to morphological cluster H5 and large sized individuals of genetic cluster 4.

Despite some evidence for phenotypic-environment associations and drift-migration equilibrium conditions, we found only weak support for positive relationships between genetic and ecological or morphological differentiation and hence a putatively adaptive basis for occupying distinct habitat niches. Although some genetic clusters varied slightly in the frequencies of different morphologies, they all contained individuals exhibiting all three body types and all five head shape types. Furthermore, colouration and depth associations were not absolute: all five clusters contained individuals with each identified main colour and were captured in areas of differing fish-community structure. When compared to collated data from other large, postglacial lakes harbouring lake trout morphs, such relatively weak or little association between morphological and genetic differentiation appears common (Table 1).

Table 1

A comparison of lake physical attributes and sympatric differentiation of lake trout as currently described in the literature.

	ML (this study)	GBL	GSL	LS	AL
Surface area (km2)	2,335[23]	31,328[20]	28,568[22]	82,100[14]	792[21]
Depth; mean/max (m)	75 / 183[75]	90 / 450[20]	73 / 614[76]	147 / 406[14]	86 / 283[77]
Secchi depth (m)	9[72]	20–30[76,78]	~ 9[76]	15–20[79]	10[80]
# of basins	2	5[20]	2[22]	2	4[21]
# of fish species	12	15[20]	21[14]	87[14]	10[80]
Shoreline distance (km)	1,509 (+ 967 isl.)	2,719 (+ 824 isl.)[18,78]	3,057[14]	2,938[14]	391 (+ 214 isl.)
# of mtDNA lineages	4[26]	1[81]	3[26]	4[26]	2[21]
# of genetically diff. clusters	5	2[19]	Not tested to date	3[17]	3‡[21]
# of morph clusters	3 body, 5 head	3–4[20,81]	3[14]	4[14]	2[21]
Ass. btw morph and genetic?	Weak	None[20]–some[19]	Not tested to date	Some[17]	None[21]
Genetic diff. (global FST)	0.017–0.028	0.008[19]	Not tested to date	0.024–0.033[17]	0.022[21]
AR¯ (# of loci)	29.51 (19)	8.30 (22)[19]	Not tested to date	3.01 (9)[14]	11.8 (8)[21]
‘Shallow’ body morphology (0-30m)	1. Lean, streamlined, dark (black, brown)	1. Streamlined, small, dark/silvery2. ‘Lean-like’, large, silvery3. Stream-lined, dark, red fins4. ‘Redfin-like’, darker, red fins[14]	1. 'Lean-like', streamlined, large, light*[22]	1. ‘Lean’: large-bodied, silvery or light[14,60]	1. ‘Lean-like’ morphs (streamlined, silvery)*[21]
‘Mid-depth’ body morphology (30-100m)	1. Large, deep-bodied, various colouration	No data to date	1. ‘Siscowet-like’, deep-body, large, dark†2. Deep body†[14,22]	1. ‘Humper’: small, narrow-bodied, silvery or light colour2. ‘Redfin’: robust, large/deep-bodied, dark, red/yellow fins[14,60]	1. ‘Siscowet-like’, deeper body§[21]
‘Deepwater’ body morphology (>100m)	1. Small, deep-bodied, light or silvery	No data to date	No LT captured	1. ‘Siscowet’: large, deep-bodied, dark colour[14,60]

Mistassini Lake (ML), Great Bear Lake (GBL), Great Slave Lake (GSL), Lake Superior (LS) and Atlin Lake (AL).

* Captured between 0–50 m

†Captured between > 50–100 m

‡ Three genetically-differentiated populations and their descriptive statistics were described in Northrup et al. [21] for the Atlin-Tagish Lake system (two lakes interconnected by river in which lake trout are also found)

§captured between 70–150 m.

Factors influencing low genetic differentiation, comparisons with other postglacial fishes

Where striking examples of multiple, sympatrically-occurring morphs of lake trout have been described, they are almost universally found in very large postglacial lakes [82]. This specificity for large lakes probably relates to (i) the large Ne that a large lake size confers towards sustaining multiple evolving populations of a large-bodied predator fish, and (ii) the different habitat niches available in larger lakes that lake trout can exploit, including shallow and deep water habitat, extensive shoreline distance, multiple basins and multiple prey fish species (Table 1). Notably, however, genetic differentiation among documented morphs is always low or relatively weak (range of global FST of 0.008–0.022; [17,19,21]; this study), and relative to previously studied fishes in Mistassini Lake (brook trout, global FST = 0.071; walleye, global FST = 0.045; [25,83]). At least four specific aspects of the ecology, behaviour and genetics of lake trout might facilitate such low differentiation.

Reduced natal fidelity and lifecycle complexity

Compared to the tributary spawning nature of brook trout and walleye, which generates considerable spatial isolation (tens to hundreds of kilometres) and a more complex lifecycle, Mistassini lake trout spawn along lake shorelines, and local fishing experts have described the spawning sites of various morphs in close geographic proximity (K. Marin 2013 pers. comm.). Although natal fidelity is well documented among postglacial lake-dwelling salmonids [84], the evidence for natal fidelity in lake trout is weak and even contradictory [85,86]; the close proximity of spawning site locations in Mistassini Lake may therefore facilitate gene flow and disfavour reproductive isolation in the long-term.

Constraints from gene flow

Consideration of Ne and m revealed significant asymmetries in m among genetic clusters, and hence the possibility that gene flow constrains adaptive divergence [87] in certain habitat niches for lake trout. Namely, m was significantly higher from the large genetic cluster occupying deeper waters (#1) to large and small clusters occupying mid-depths and shallow waters (#2, 4, 5). Adaptive divergence for occupying some habitat niches in shallower waters may thus be constrained; alternatively, selection against migrants may be stronger against individuals from shallower than deeper genetic clusters.

Constraints from phenotypic plasticity

Whereas defined genetic clusters exhibited summer depth preferences, individuals from different clusters were captured throughout the water column. The striking morphological variation was more continuous within than among genetic clusters, suggesting a highly plastic propensity for occupying multiple habitat niches. Habitat use among individual lake trout has indeed been found to be highly variable but consistent between years [88-90].

Phenotypic plasticity is thought to both facilitate and constrain adaptive divergence and reproductive isolation during the process of ecological speciation [2,91]. Sympatric differentiation in lake trout might be quite unique in exemplifying how plasticity may constrain further steps towards ecological speciation. For example, in Mistassini Lake, low genetic differentiation was unrelated to a potential downward bias in FST and a lack of migration-drift equilibrium, but rather, related to a high contemporary gene flow and large Ne. And, in contrast to other cases in north temperate fishes [92,93], including co-generic Arctic charr [8], remarkable phenotypic variation in Mistassini lake trout and in other lakes does not result in appreciable genetic segregation nor a clear correspondence between different phenotypes and different habitat niches, despite evidence for the availability of the latter. Such discontinuous adaptive variation with minor reproductive isolation implies that i) selection against migrants or hybrids between diversifying lake trout populations is not very strong, ii) morphological plasticity is not very costly in this species, and/or iii) the extent of adaptive divergence might be constrained by adaptive plasticity [2].

If low genetic differentiation in sympatric lake trout is constrained by plasticity, then a key question is why plasticity might more readily constrain later steps of ecological speciation in this species relative to congeneric Arctic charr. One possibility is that the evolution of lake trout was initially born out of specialization for terminal piscivory in cold-water lakes [82], whereas Arctic charr do not exhibit such near ubiquity in trophic specialization [6]. Such specialization might act as a developmental limitation to further diversification through phylogenetic niche conservatism, and may in part explain the rarity of genetically-distinguishable morphs in lake trout. A second possibility is that species differences in life history might have evolved in response to previous glacial cycles. The evolution of lake trout, unlike that of Arctic charr, was refined to the freshwater environment as opposed to marine and freshwaters [82]. The evolution of the anadromous life history stage of Arctic charr may have therefore played a critical role in the high precision of natal homing in this species [94]. Hence, adaptive plasticity and adaptive divergence might be more coupled, increasing the likelihood of greater genetic segregation arising where different habitat niches exist. Collectively, the unique evolutionary history of two congeneric species may allow plasticity to act as a constraint to ecological speciation at early stages in one species (lake trout) but help facilitate it at later stages in another (Arctic charr).

Genetic diversity and evolutionary origins of sympatric lake trout morphs

This study on the complete spatial distribution of Mistassini lake trout (n = 636) and an independent study that screened Mistassini lake trout of unknown location (n = 46) using the same 19 microsatellite loci [95], have found correspondingly high genetic diversity (mean NA = 29.51 and 17.26, respectively). Indeed, we found no evidence that Mistassini lake trout have experienced bottlenecks (S9 File) and the pristine nature of the lake probably means that the large Ne of inferred genetic clusters has not been reduced since colonization. Mistassini morphs therefore display high genetic diversity relative to morphs in other large, postglacial lakes (mean AR = 3.01–11.80; Table 1). However, studies on other lakes used fewer (exception: Great Bear) and different microsatellite loci. Furthermore, lower genetic diversity in other lakes may be due to incomplete spatial sampling, lower numbers of colonizing mtDNA lineages (one in Great Bear vs. up to four in Mistassini [19,26]; or historical reductions from overfishing, invasive species and/or intensive stocking (i.e., Lake Superior) [96].

We found that genetic differentiation among clusters accounting for the mutational properties of microsatellites (RST) exceeded global FST (S3 File), and hence affirmed previous work [26] that multiple lake trout lineages colonized Mistassini Lake. Indeed, Mistassini’s age is too young (8,000 years) and the number of generations passed are too few (667; 12 year generation time; [28,29]) to have generated this mutationally-driven differentiation entirely within sympatry. Yet an equilibrium between drift and gene flow has likely been reached since lake colonization, due to relatively high contemporary gene flow compared to Ne (Fig S7.1; Table S7.2 in S7 File). Because at least one lake (Great Bear) was postglacially colonized by a single mtDNA lineage [19,26], an allopatric phase does not appear to be a prerequisite for sympatric differentiation in lake trout as observed in some postglacial fishes (e.g., [97,98]). However, as our study reveals, individuals from different lake trout lineages have a propensity to occupy and exploit multiple available habitat niches.

Conclusion: early-stage ecological speciation and conservation implications

In several large, postglacial lakes, lake trout exhibit non-random structuring or complexes of weakly diversifying populations with often striking phenotypic variation associated or partially associated with depth, spatial distribution and/or resource exploitation (Fig 5; Table 1); these morphs do not easily correspond to the classical benthic-limnetic morphs found in other postglacial fishes (e.g., [2,5,97]). In a remarkably pristine system (Mistassini Lake), our study revealed likely factors influencing the progress or lack of progress at earlier stages of population differentiation in the context of ecological speciation: genetic clusters appear to harbour considerable phenotypic plasticity, their habitat use distinctions were not clearly defined with individuals patchily occupying different habitat niches, reproductive habitats of different morphs were in close spatial proximity, and any adaptive divergence generated from differential habitat use appears to only generate weak reproductive isolation. Plasticity in particular, may act as a constraint to further diversification and place lake trout at the early-stages of the ecological speciation continuum [2].

Fig 5

Morphological variation among individual lake trout and across all genetically-demarcated clusters within Mistassini Lake, Québec.

Coloured circles and numbers represent the genetically-demarcated clusters that individual lake trout were assigned to.

Whether lake trout morphs have ‘plateaued’ in their steps towards true ecological speciation, or whether further adaptive divergence and strengthening of reproductive isolation is still to come in a future context, the maintenance of such differentiation may be important for species persistence in the face of changing environments [5]. Indeed, lake trout are found in tens of thousands of lakes but sympatric morphs have been documented in only a handful. Our research suggests that different morphs may be more or less susceptible to different harvesting techniques (gillnetting vs. angling). Their occupancy of large, deep lakes means that such lakes act as reservoirs of genetic diversity and cold water refuge habitat for the species in the face of future climate change. As a terminal predator, conservation plans for lake trout or for large postglacial ecosystems may benefit from a joint consideration of the role and survival of both prey, competitors (e.g., brook trout, walleye, whitefish, etc.) and conspecifics that thrive across a diversity of habitats, to maintain the processes generating biologically significant differentiation.

Loci and polymerase chain reaction (PCR) information.

(PDF)

Click here for additional data file.

Supplemental genetic clustering approaches.

(PDF)

Click here for additional data file.

Genetic differentiation by stepwise mutation (RST).

(PDF)

Click here for additional data file.

Supplemental morphological figures and tables.

(PDF)

Click here for additional data file.

Multinomial regression analysis.

(PDF)

Click here for additional data file.

Supplemental genetic data of lake trout in Mistassini Lake.

(PDF)

Click here for additional data file.

Effective population size (N) and gene flow Mistassini lake trout.

(PDF)

Click here for additional data file.

Supplemental tables and figures for associations between genetic, morphological and ecological differentiation.

(PDF)

Click here for additional data file.

Assessment of recent bottlenecks of Mistassini lake trout.

(PDF)

Click here for additional data file.