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).