Results:

Ecomorph divergence in head shape:

For head shape the benthivore and planktivore ecomorphs were separated across PC1 (31% variance explained), except from Loch Dughaill where the ecomorphs separate along PC2 (17.3%) (Figure 1C). Individuals with a positive PC1 score had shallower heads with larger eyes than those with negative PC1 scores (Figure 1E). For PC2, a more positive score suggested a longer head shape. The benthivore ecomorphs generally have a more negative PC2 score suggesting their heads are shorter than the planktivore ecomorphs (Figure S2).
We compared the magnitude and direction of phenotypic change for head shape between the ecomorphs across pairs, to determine how similar the divergences were, through a phenotypic trajectory analysis (PTA). We found that for all pairwise comparisons, the angle of difference in phenotypic trajectories (θ) was significantly different (P < 0.05) (Table 1). Similarly, almost all differences in trajectory lengths between ecomorphs pairs (ΔL) were significantly different with the exception of the Awe vs na Sealga comparison. These results suggest that head shape morphology is variable across lakes. This is in agreement with our ANOVA model (PC ~ ecomorph + lake + ecomorph x lake) which found that the ecomorph x lake interaction term explained most variation (η2Eco*Lake = 0.565) suggesting that the effect of lake environment and/or evolutionary history strongly impacts the direction and magnitude of head shape divergence between ecomorphs, and that head shape is not strictly parallel across lakes (Table S3).

Genomic regions associated with head shape:

To determine the genomic variation underpinning head shape, we performed a genome-wide association analysis (GWAS) on a set of 13,071 SNPs (Figure S3). Using the PC scores from the head shape analysis (PCs 1-5), a Bayesian Sparse Linear Mixed Model (BSLMM) showed that the proportion of phenotypic variance in head shape explained by genetic variation (PVEHead) in the SNP dataset was 0.62 with the proportion of phenotypic variance explained by large (“sparse”) effect loci (PGEHead) was 0.82. This is supported by the ⍴ (rho) value for head shape that suggests that the head shape phenotype is controlled largely by a few large-effect loci (⍴Head= 0.792).
Applying Linear Mixed Models to identify SNPs highly associated with head shape variation found a total of 82 SNPs (66 SNPs mapped on 27 of 39 chromosomes, 16 mapped to unanchored scaffolds) that showed a significant association with variation in head shape (Bonferroni corrected P-value < 0.05; Fig. 2A) with these SNPs broadly distributed across the genome.

Genomic differentiation at SNPs associated with head shape:

We investigated whether these head shape-associated SNPs were highly diverged between the ecomorphs in all lakes, consistent with a shared genomic bases for these phenotypes, or whether they were specific to certain populations suggesting the deployment of different genetic pathways leading to the similar phenotypes across pairs. We found a total of three SNPs that were diverged in all four lakes (Fig. 3A).
We aimed to identify if those SNPs associated with head shape showed signs of response to divergent or positive selection in all four lakes. Mean genetic differentiation (FST) and absolute divergence (DXY) between ecomorphs in lochs Dughaill and Tay were elevated amongst associated SNPs when compared to the background (Figure 4, 5) (Table S4). FST and DXY were significantly lower than background between ecomorphs in Loch Awe for the associated SNPs. There was no significant difference in FST or DXY between associated SNPs and the background in na Sealga. These results suggest that the SNPs associated with head shape are not similarly responding to, or are under selection, across all lakes. These patterns were not influenced by linkage, because recombination rates for genomic regions around head shape-associated SNPs did not differ significantly from genomic background (Figure S4).

Genes and QTLs associated with head shape variation:

Focusing on the location of the head shape-associated SNPs relative to known genes in the charr genome, we found that 38 of the 82 SNPs were located within annotated genes (Table S5). GO term analyses found that the genes containing head-shape associated SNPs showed over-enrichment for GO terms related to odontogenesis (GO:0042476), cranial skeleton system development (GO:1904888) among other processes and functions (Table S6) when compared to all genes containing SNPs in our dataset.
To examine if any of these genomic associations were shared across other species, we compared the positions of the head shape-associated SNPs to QTL markers from across salmonid species. We found that three of the head shape-associated SNPs were found within ±100,000 bp of the peak positions of two mapped QTLs (Table 2). These two QTLs were previously found to be associated with body shape morphology in lake trout and lake whitefish .

Ecomorph divergence in body shape:

For body shape, all four ecomorph pairs showed separation along PC1 (30.3%) (Figure 1D) however, the pair from Loch Dughaill diverged in a different direction, along PC2 (24.9%). Individuals with a positive PC1 score (e.g., the benthivore morphs at Awe, na Sealga, and Tay) have shallower, more elongated body shapes (Figure 1F). The patterns across PC2 suggest that a more positive score is associated with a deeper body (Figure S5).
In the PTA, we again found that all differences in the magnitude of phenotypic change between the pairs (ΔL) was significant with the exception of the Awe-na Sealga comparison (Table 1). When comparing angle of difference in trajectories (θ), we found all angles between significant with the exception of the na Sealga-Tay comparison suggesting body shape may show some parallelism across lakes. The ecomorph term in the ANOVA model (η2Eco = 0.301) for body shape explained more variation than the ecomorph x lake term (η2Eco*Lake = 0.135), suggesting the ecological niche is a stronger determinant of body shape, indicating some level of phenotypic parallelism across lakes. However, the unique evolutionary history (the lake term) explained most variation (η2Lake = 0.414) suggesting that the absolute position of ecomorph pairs in the multivariate space differs between lakes and is strongly influenced by differences in the evolutionary histories or environment of each of the pairs (Table S3).

Genomic regions associated with body shape:

A BSLMM found that the proportion of phenotypic variance in body shape explained by the SNP dataset (PVEBody) was 0.82 and the proportion of phenotypic variance explained by large (“sparse”) effect loci (PGEBody) was 0.39 and the ⍴ (rho) value (⍴Body) was 0.425. By analysis with LMMs, we found 180 SNPs significantly associated with body shape variation (144 SNPs mapped to 34 chromosomes, 36 SNPs mapped to unplaced scaffolds) (Figure 2B). We found that 10 of these SNPs were present in all four ecomorph pairs (Figure 3B) and broadly distributed across the genome (on 7 chromosomes).

Genomic differentiation at SNPs associated with body shape:

For FST, we found that the body shape-associated SNPs had a higher mean value than the background at Loch Dughaill and Tay (Figure 4, 5, Table S4) but no notable difference at Loch Awe or na Sealga. DXY was significantly higher at Tay for the associated SNPs while it was significantly lower at Loch Awe. There was no significant difference in DXY between associated SNPs and the background at Loch Dughaill and naSealga. The mean recombination rate in regions containing body shape SNPs did not differ from the background (Figure S4).

Genes and QTLs associated with body shape:

Relative to known genes in the Salvelinus sp. genome, 89 of the body shape-associated SNPs were located within annotated genes (Table S5). The body shape-associated genes showed overrepresentation for genes involved in skeletal system (GO:0001501), face (GO:0060324), eye (GO:0060041, GO:0001745), and mouth development (GO:0060021) among other processes and functions (Table S6). We found five SNPs in close proximity to four known QTLs (Table 2). These QTLs were previously found to be associated with body shape morphology in lake whitefish, body weight in Arctic charr and body shape morphology in lake trout .

Comparisons between head shape and body shape:

We found in our PTA that comparisons in head shape showed greater mean differences in the magnitude and direction of phenotypic change (ΔLHead = 0.053 ± 0.035 s.d., mean θHead= 69.87° ± 17.54 s.d. mean) compared to body shape (mean ΔLBody = 0.016 ± 0.011 s.d., mean θBody= 61.62° ± 22.63 s.d.) (Table S7). Both our PTA and ANOVA suggest that body shape shows slightly stronger patterns of phenotypic parallelism than head shape, however, both phenotypes show substantial deviations from strict parallelism across lakes.
Our association analyses indicated a significant shared genetic basis behind body shape and head shape. 50 of the SNPs found in our study appeared associated both with head and body shape morphology (212 SNPs identified: 32 SNPs associated with head shape, 130 associated with body shape, and 50 associated with head and body shape), which exceeds random expectation (hypergeometric test; P = 6.633e-16). Of these head- and body-shape shared SNPs, 38 mapped to 20 chromosomes and 12 mapped to unplaced scaffolds (Figure 2). These SNPs show overrepresentation for terms related to brain (GO:0030900, GO:0021575) and heart development (GO:0003007), and regulation of cell shape (GO:0008360) among other processes (Table S6). With a number of SNPs shared between both head and body shape, two of the QTLs we identified as near associated SNPs were near SNPs shared for both phenotypes. These QTLs related to body shape in lake trout and lake whitefish respectively (Table 2) .
This shared genetic basis is reflected in the PTA, which showed that there was a positive linear relationship when comparing trajectory lengths for head shape and body shape across lakes. Specifically, pairs in which the ecomorphs have diverged to a similar extent in head shape have also diverged to a similar extent in body shape (adjusted R2 = 0.99, P < 0.001) (Figure 6). A similar positive relationship was seen when comparing differences in the direction of phenotypic change, although this was non-significant (adjusted R2 = 0.29, P = 0.154) (Figure S6).