\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,greek,polish,english]{babel}



\usepackage{float}






% Edit this header.tex file to include frontmatter definitions and global macros

% Add here any LaTeX packages you would like to load in all document blocks
% \usepackage{xspace}

% Add here any LaTeX macros you would like to load in all document blocks
% \def\example{This is an example macro.}

% -----

\iflatexml
% Add here any LaTeXML-specific commands

% -----

\else
% Add here any export style-specific LaTeX commands. These will only be loaded upon document export.
% \paperfield{Subject domain of my document}
% \keywords{keyword1, keyword2}
% \corraddress{Author One PhD, Department, Institution, City, State or Province, Postal Code, Country}
% \fundinginfo{Funder One, Funder One Department, Grant/Award Number: 123456.}
\fi


\begin{document}

\title{Genomic underpinnings of head and body shape in Arctic charr ecomorph
pairs}



\author[1]{Sam Fenton}%
\author[1]{Arne Jacobs}%
\author[2]{Colin Bean}%
\author[3]{Colin Adams}%
\author[3]{Kathryn Elmer}%
\affil[1]{University of Glasgow College of Medical Veterinary and Life Sciences}%
\affil[2]{Affiliation not available}%
\affil[3]{University of Glasgow}%


\vspace{-1em}



  
  \date{March 13, 2025}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Across its Holarctic range, Arctic charr (Salvelinus alpinus)
populations have diverged into distinct trophic specialists across
independent replicate lakes. The major aspect of divergence between
ecomorphs is in head shape and body shape, which are ecomorphological
traits reflecting niche use. However, whether the genomic underpinnings
of these parallel divergences are consistent across replicates was
unknown but key for resolving the substrate of parallel evolution. We
investigated the genomic basis of head shape and body shape morphology
across four benthivore-planktivore ecomorph pairs of Arctic charr in
Scotland. Through genome-wide association analyses, we found genomic
regions associated with head shape (89 SNPs) or body shape (180 SNPs)
separately and 50 of these SNPs were strongly associated with both body
and head shape morphology. For each trait separately, only a small
number of SNPs were shared across all ecomorph pairs (3 SNPs for head
shape and 10 SNPs for body shape). Signs of selection on the associated
genomic regions varied across pairs, consistent with evolutionary
demography differing considerably across lakes. Using a comprehensive
database of salmonid QTLs newly augmented and mapped to a charr genome,
we found several of the head and body shape associated SNPs were within
or near morphology QTLs from other salmonid species, reflecting a shared
genetic basis for these phenotypes across species. Overall, our results
demonstrate how parallel ecotype divergences can have both
population-specific and deeply shared genomic underpinnings across
replicates, influenced by differences in their environments and
demographic histories.%
\end{abstract}%



\sloppy


\section*{Genomic underpinnings of head and body shape in Arctic charr
ecomorph
pairs}

{\label{genomic-underpinnings-of-head-and-body-shape-in-arctic-charr-ecomorph-pairs}}

Sam Fenton\textsuperscript{1}

https://orcid.org/0000-0001-5011-0329

Arne Jacobs\textsuperscript{1}

http://orcid.org/0000-0001-7635-5447

Colin W. Bean\textsuperscript{1,2}

https://orcid.org/0000-0003-3502-0995

Colin E. Adams\textsuperscript{1,3,}*

https://orcid.org/0000-0003-2470-9754

Kathryn R. Elmer\textsuperscript{1,4,}*

https://orcid.org/0000-0002-9219-7001

1 School of Biodiversity, One Health \& Veterinary Medicine, University
of Glasgow, Glasgow, UK.

2 NatureScot, Clydebank, UK.

3 Scottish Centre for Ecology and the Natural Environment, University of
Glasgow, Glasgow, UK.

4 corresponding author: Kathryn.Elmer@glasgow.ac.uk

* joint senior authors

\section*{Running title: Genomics of charr head and body
shape}

{\label{running-title-genomics-of-charr-head-and-body-shape}}

\section*{Abstract:}

{\label{475102}}

Across its Holarctic range, Arctic charr (\emph{Salvelinusalpinus} )
populations have diverged into distinct trophic specialists across
independent replicate lakes. The major aspect of divergence between
ecomorphs is in head shape and body shape, which are ecomorphological
traits reflecting niche use. However, whether the genomic underpinnings
of these parallel divergences are consistent across replicates was
unknown but key for resolving the substrate of parallel evolution. We
investigated the genomic basis of head shape and body shape morphology
across four benthivore-planktivore ecomorph pairs of Arctic charr in
Scotland. Through genome-wide association analyses, we found genomic
regions associated with head shape (89 SNPs) or body shape (180 SNPs)
separately and 50 of these SNPs were strongly associated with both body
and head shape morphology. For each trait separately, only a small
number of SNPs were shared across all ecomorph pairs (3 SNPs for head
shape and 10 SNPs for body shape). Signs of selection on the associated
genomic regions varied across pairs, consistent with evolutionary
demography differing considerably across lakes. Using a comprehensive
database of salmonid QTLs newly augmented and mapped to a charr genome,
we found several of the head and body shape associated SNPs were within
or near morphology QTLs from other salmonid species, reflecting a shared
genetic basis for these phenotypes across species. Overall, our results
demonstrate how parallel ecotype divergences can have both
population-specific and deeply shared genomic underpinnings across
replicates, influenced by differences in their environments and
demographic histories.

Keywords: parallel evolution, ecomorphology, QTL, genomics, fish,
population genomics

\section*{Introduction:}

{\label{introduction}}

Parallel evolution is an evolutionary process and outcome by which
similar phenotypes arise and establish in multiple independent
populations in separate environments (Bolnick et al., 2018; Elmer \&
Meyer, 2011; Schluter, 1996). These replicate evolutionary outcomes
suggest that similar environments impose similar selective pressures on
organismal phenotypes, with only a small number of phenotypic solutions
favoured in that context. While the existence of parallel phenotypes is
well established in natural populations (Oke et al., 2017; Siwertsson et
al., 2013), the extent to which those are associated with similarly
shared genomic underpinnings is rarely examined (Magalhaes et al., 2021;
McGirr \& Martin, 2018; Ravinet et al., 2016). Indeed, the appearance of
the same phenotypes through parallel evolution processes does not mean
that the same genomic processes underpin those similar evolutionary
outcomes across replicates (Conte et al., 2012; Elmer \& Meyer, 2011).
Similar phenotypic outcomes could result from alternative genetic
pathways. This might arise because of differing demographic histories
(Elmer et al., 2014), variable genetic backgrounds (Arendt \& Reznick,
2008; Kowalko et al., 2013) or the involvement of alternative splicing,
differential gene expression or post-translational modifications
resulting in phenotypic parallelism (Filteau et al., 2013; Jacobs \&
Elmer, 2021; McGirr \& Martin, 2018).

One classic example of parallel evolution is the replicated divergence
across northern freshwater lakes of distinct trophic specialists, also
known as ecomorphs or ecotypes. These occur abundantly in salmonid
fishes in recently glaciated lakes, such as lake whitefish
(\emph{Coregonus clupeaformis}) (Landry et al., 2007), lake trout or
lake charr (\emph{Salvelinus namaycush}) (Baillie et al., 2016), and
Arctic charr (\emph{Salvelinus alpinus}) (Doenz et al., 2019; Jensen et
al., 2017). In Arctic charr, ecomorphs associated with divergence along
the depth axis and ecological niche, typically forming pelagic and
benthic foraging specialists (Elmer, 2016; Klemetsen, 2010). Two key
traits involved in this divergence are head shape and body shape, both
having functional significance; head shape being important for foraging
and prey specialisation and body shape important for swimming behaviour
and niche use (Adams \& Huntingford, 2002; Skoglund et al., 2015).
Additionally, ecomorphs often differ in other complex traits such as
body size, colouration and spawning time (Gardu\selectlanguage{ngerman}ño-Paz et al., 2012;
Jonsson \& Jonsson, 2001).

Previous research on benthivorous-planktivorous ecomorph pairs of Arctic
charr has shown that the repeatability into ecomorph is highly parallel
and that individual traits related to head shape also show some
parallelism (Adams et al., 2008; Jacobs et al., 2020). Genomic analysis
to-date suggested limited genetic parallelism across lakes, with many
differences in demography and colonisation history across the breadth of
charr distribution (Brachmann et al., 2022; Jacobs et al., 2020;
Salisbury et al., 2020). However, this previous research on parallelism
between Arctic charr ecomorphs focused on the patterns of the genomic
response to selection such as outlier loci, which are known to have high
rates of false-positives and false-negatives and are indirect approaches
to understanding morphology, strongly influenced by evolutionary genomic
and demographic histories (Ravinet et al., 2017). Therefore, we lack
knowledge of the genomic regions associated with key phenotypic
differences in charr and the extent to which these are shared across
ecomorph pairs (Elmer, 2016).

Because of their vital role in foraging and swimming, the genetics of
head and body morphology has been explored as QTL studies in many fish
species, including salmonids (Boulding et al., 2008; Küttner et al.,
2014; Laporte et al., 2015; Smith et al., 2020). For example, in lake
whitefish, as many as 138 different quantitative trait loci (QTLs)
related to body shape have been identified, implying that the trait is
highly polygenic in this species (Laporte et al., 2015). In other
species, body shape has been found to be less polygenic, with a small
number of QTLs identified related to benthic-limnetic differences
between ecomorphs of lake trout (Smith et al., 2020) and a single region
that differentiates ecomorphs in Atlantic cod (\emph{Gadus morhua)}
(Hemmer-Hansen et al., 2013). However, the extent to which parallel
phenotypes have the same genomic underpinnings across replicates seems
to vary greatly between species and locality. Studies in three-spine
stickleback (\emph{Gasterosteus aculeatus}), rainbow trout
(\emph{Oncorhynchus mykiss}), and sockeye salmon (\emph{Oncorhynchus
nerka}) all suggest that while some genetic variation that underlies
morphological variation can be shared across replicates, this is rarely
the case across the species' entire range and the degree of genetic
parallelism that is shared is often low (Fang et al., 2020; Larson et
al., 2019; Weinstein et al., 2019). Often only a few key genes or loci
are shared across replicates, and in some cases, they underlie morph
differentiation in multiple species (Jacobs et al., 2017; Salisbury \&
Ruzzante, 2021).

In this study, we investigated the genomic underpinnings of head and
body shape morphology in replicate ecomorph pairs of Arctic charr across
four independent lakes in Scotland. First, we determined the extent of
phenotypic parallelism in head and body shape morphology across ecomorph
pairs. We examined these two traits separately as they are known to
involve different axes of ecological specialisation and can have
different genetic bases (Boulding et al., 2008; Küttner et al., 2014;
Smith et al., 2020). Second, we investigated the genome-wide
underpinnings of these phenotypes by identifying SNPs strongly
associated with head shape variation and body shape variation. We then
evaluated the genomic organisation of these head shape- and body
shape-associated SNPs, specifically to determine if they were
distributed widely across the genome, co-localised in genes or genomic
regions, or within known salmonid QTLs. To do this, we updated and
augmented an existing QTL database (Jacobs et al., 2017) and mapped this
to the orthologous location in the \emph{Salvelinus sp.} genome. Third,
we examined evolutionary genomic background by analysis of selection on
those loci we found to be associated with head shape and with body
shape. Shared signals of selection and elevated differentiation and
divergence in independent ecomorph pairs would suggest similar processes
across evolutionary replicates. Finally, we synthesised the findings to
identify shared genomic underpinnings and similarities in phenotypic
divergence between head and body shape, as two interlinked components of
ecomorph.

\section*{Materials and Methods:}

{\label{materials-and-methods}}

\subsection*{Study populations:}

{\label{study-populations}}

Fish populations were examined from four lakes in Scotland: Loch Awe,
Loch Dughaill, Loch na Sealga, and Loch Tay (Figure 1A). Each lake
contains two ecomorph populations, one benthivorous and the other
planktivorous, with each lake representing an ecomorph pair (Garduño-Paz
et al., 2012; Hooker et al., 2016; Jacobs et al., 2020). Explictily, we
investigated the simiarlities in divergences of head and body shape
between the two different ecomorphs across pairs. Short-read genomic
data (ddRADseq NCBI short read bioproject: PRJNA607173) and fish
photographs were drawn from previous research (Jacobs et al., 2020) and
reanalysed here.

\subsection*{Morphological analyses:}

{\label{morphological-analyses}}

To investigate patterns in head morphology and body shape separately,
landmarks for head (16 landmarks) and body shape (14 landmarks) were
placed using \emph{ImageJ v1.50i} (Schneider et al., 2012) \emph{}
(Figure 1E-F). To allow us to cover the whole-body shape of each
individual, several markers related to head shape were retained in the
body shape analysis. As such, the head and body shape analyses are not
completely independent but are focused on each specific phenotype and
use a distinct set of landmarks. In total, 341 individual fish were
landmarked for head shape (Awe\_Bn=37, Awe\_Pl=31, Dughaill\_Bn=43,
Dughaill\_Pl=54, naSealga\_Bn=42, naSealga\_Pl=20, Tay\_Bn=41,
Tay\_Pl=73) and 335 were landmarked for body shape (Awe\_Bn=34,
Awe\_Pl=30, Dughaill\_Bn=43, Dughaill\_Pl=54, naSealga\_Bn=42,
naSealga\_Pl=19, Tay\_Bn=40, Tay\_Pl=73).

General Procrustes Analysis (GPA) was performed to standardise each
individual's shape for size and orientation using \emph{geomorph v3.0.7}
(R package) (D. C. Adams \& Otárola-Castillo, 2013). Following this
procedure, shape data were standardised for any residual size effects
using the log of centroid size to correct for allometry. Size corrected
data then were used in all subsequent analyses. Principal Component
Analyses (PCA) were conducted on head shape and body shape separately
using the \emph{plotTangentSpace} and \emph{arrayspecs} functions and
plots were made in \emph{ggplot2} \emph{v3.3.5} (R package) (R Core
Team, 2021; Wickham, 2016).

Analysis of variance (ANOVA) models were used to determine the relative
contributions of parallel and non-parallel aspects of the morphological
divergence across ecomorph pairs (Langerhans \& DeWitt, 2004). Rather
than using a single PC (Principal Component) score for the ANOVAs (ANOVA
model: PC \textasciitilde{} ecomorph + lake + ecomorph x lake), we used
a combined weighted PC variable with scores from multiple PCs to account
for a higher proportion of biologically important variance in a single
variable. For this variable, we combined all consecutive PCs that
explained more than 5\% of the total variance in shape. For both head
shape and body shape, this was PC1 to PC5 (Figure S1). To weight the
scores within the combined variable, each PC score was multiplied by the
amount of variance explained by that PC (i.e., each individual PC1 score
x proportion of variance explained by PC1) before all five weighted
scores were summed for each individual. The \emph{EtaSq} function in the
\emph{BaylorEdPsych v0.5} (R package) (Beaujean, 2012) was used to
estimate the effect size of each model term.~In our model, the ecomorph
term represents the parallel or shared term, the ecomorph*lake
interaction is the non-parallel (non-shared) term, and the lake term
represents unique evolutionary history. ANOVA results for PC1 through
PC5 separately can be found in Table S3.

Phenotypic Trajectory Analyses (PTA) were performed on the procrustes
scores using the \emph{trajectory.analysis} function in \emph{geomorph}
to look at the extent and direction of phenotypic change between
ecomorphs. Magnitude of divergence is described by the length of
trajectories (L) while the angle between trajectories (\selectlanguage{greek}θ\selectlanguage{english}) describes
their direction in phenotypic space.~ This approach allows us to
determine how parallel the trajectories of each ecomorph pair are to one
another by using the difference in phenotypic trajectory length
(\selectlanguage{greek}Δ\selectlanguage{english}L\textsubscript{P}) and the direction of phenotypic trajectories
(\selectlanguage{greek}θ\selectlanguage{english}\textsubscript{P}). The significance of differences in trajectory
lengths and differences in trajectory direction were assessed using
1,000 permutations (Adams \& Collyer, 2013; Stuart et al., 2017).
Phenotypic divergence between ecomorphs in different lakes was
considered to be parallel if the direction and/or magnitude of phenotype
change did not differ significantly between the pairs (P \textless{}
0.05).

\subsection*{Population genomics:}

{\label{population-genomics}}

Filtered 75bp reads for each individual, generated via ddRADseq from
Jacobs \emph{et al}. (2020) and accessed from (ddRADseq NCBI short read
bioproject: PRJNA607173), were mapped using \emph{bwa mem} and
\emph{SAMtools} using settings described in that paper to the
\emph{Salvelinus} sp\emph{.} genome from NCBI (ASM291031v2). The number
of reads per individual ranged from 1 to 3.5 million. RAD loci were
built in the \emph{gstacks} module of \emph{Stacks v2.53} (Rochette et
al., 2019) for 200 individuals (Awe\_Bn= 26, Awe\_Pl=29,
Dughaill\_Bn=28, Dughaill\_Pl=27, naSealga\_Bn=18, naSealga\_Pl=20,
Tay\_Bn=21, Tay\_Pl=31). SNPs were retained in the \emph{populations}
module of \emph{Stacks} if they met the following criteria: present in
66\% of all individuals in each population and across all populations, a
minimum minor allele frequency of 0.05, maximum observed heterozygosity
of 0.5. Each ecomorph within a lake was considered to be a discrete
population. The script filter\_hwe\_by\_pop.pl to filter out sites
outside Hardy-Weinberg equilibrium within populations (available at
\url{https://github.com/jpuritz/dDocent}). \emph{vcftools v0.1.13}
(Danecek et al., 2011) was used to filter to a minimum coverage of 5x
and a maximum of 50x. A principal component analysis was performed to
identify the major axes of genetic variation using \emph{SNPRelate
v1.22.0} (R package) (Zheng et al., 2012).

\subsection*{Genotype-phenotype association
analyses:}

{\label{genotype-phenotype-association-analyses}}

To determine the association between genotypes and phenotypic variation
in head or body shape, we ran Linear Mixed Models (LMM) in \emph{Gemma
v0.98.1 (Zhou \& Stephens, 2012)}. Univariate and Multivariate LMMs with
Wald's test were run using PCs 1-5 for head shape and body shape, the
SNP dataset generated for the population genomics analyses, and lake of
origin as a co-variate (Zhou \& Stephens, 2014). Missing genotypes were
imputed using \emph{LinkImpute v1.1.4} (Money et al., 2015) and a
relatedness matrix was generated using \emph{Gemma} before running the
models. We determined significant associations using
bonferroni-corrected P-values (0.05/7329 unlinked SNPs) from the Wald's
tests. The number of unlinked SNPs was determined by LD-pruning the full
SNP dataset using the snpgdsLDpruning function in \emph{SNPRelate}.
Bayesian Sparse Linear Mixed Models (BSLMM) (Zhou et al., 2013) were run
using PC1-5 variables to determine how much of the phenotypic variation
is explained by the SNPs in our dataset (PVE) and secondly how much of
that variation is explained by large-effect loci (PGE), and finally how
polygenic each phenotype is ([?]). Manhattan plots were made using
\emph{CMplot v4.0} (R package)
(\url{https://github.com/YinLiLin/CMplot}).

We subsequently determined if SNPs showing significant associations with
head shape or body shape morphology were found within annotated genes in
the \emph{Salvelinus} sp\emph{.} reference genome using \emph{BEDtools
v2.27.1 (Quinlan \& Hall, 2010)}. The functions of genes containing, or
\selectlanguage{ngerman}± 1kbp of, associated SNPs were investigated using GO term
overrepresentation analysis (ORA) and gene set enrichment analysis
(GSEA). These analyses were run using \emph{topGO v2.40.0} (R package)
(Alexa \& Rahnenfuhrer, 2020) with all genes containing any RAD loci as
the full comparison dataset. Results were summarised using \emph{REVIGO}
(Supek et al., 2011) before visualisation in \emph{Cytoscape v3.91}
(Shannon et al., 2003).

\subsection*{Comparisons to known QTLs:}

{\label{comparisons-to-known-qtls}}

Using existing information on the genetics of important phenotypes from
other salmonid species, we mapped a database of 1,338 QTL markers to the
\emph{Salvelinus} sp\emph{.} genome. This was based on a previously
published database of QTLs involved in traits related to morphology and
life history, derived from a range of salmonid species and previously
mapped to the \emph{Salmo salar} genome (Jacobs et al., 2017).
Additionally, a literature search was conducted up to April 2021 to
augment the existing database with more recently published QTLs. This
literature search was conducted in Web of Science and Google Scholar
using the search terms ``QTL'', ``quantitative trait loci'',
``salmonid'', and the common and scientific names for rainbow trout,
Atlantic salmon, Arctic charr, lake whitefish, Chinook salmon, coho
salmon, brook trout, and lake trout. QTL marker sequences were gathered
for 17 different phenotypes: body length, body shape, body weight,
Fulton's condition factor, directional change, disease resistance,
embryonic development, gill rakers, growth rate, hatching time, head
shape, parasite resistance, salinity tolerance, sexual maturation,
smolting, spawning time, upper temperature tolerance (Table S1).

Following Jacobs \emph{et al.} (2017), the strategy of mapping the
QTL-linked markers to the \emph{Salvelinus} sp. genome depended on the
QTL marker type: RAD loci were mapped using \emph{Bowtie2 v2.4.4} and
the \emph{very sensitive} pre-set; microsatellite primer sequences,
which are shorter, were mapped using \emph{Bowtie v1.3.1} allowing for 3
mismatches. QTLs for which the flanking markers mapped to different
chromosomes were removed. Redundant QTLs, i.e., where two QTLs for the
same trait from the same species mapped to the same location, were
removed only keeping the QTL with the higher PVE or LOD score (following
Jacobs et al., 2017). For QTLs where more than one marker was reported,
we attempted to map all markers. Position values for the QTLs markers
were then compared to positions of the phenotype associated SNPs using
\emph{BEDtools}, with a cut-off of ±100kbp. This value was used so that
we could consider SNPs within the range to be proximal to a QTL peak
while also accounting for the large size of many of the QTLs in the
database. In total, we successfully mapped 669 QTL-linked sets of
markers to the \emph{Salvelinus} sp. genome after removing redundant
QTLs (Table S2).

\subsection*{Genomic response to
selection:}

{\label{genomic-response-to-selection}}

We investigated if the phenotype-associated SNPs identified in our
analyses showed signals of a response to selection and if those signals
were replicated across ecomorph pairs. To test this, for each ecomorph
pair we compared F\textsubscript{ST} and D\textsubscript{XY} values for
phenotype-associated SNPs to a random background subset of SNPs. This
random subset was 100 SNPs randomly selected from the whole dataset and
the mean F\textsubscript{ST} and D\textsubscript{XY} values for those
SNPs were calculated. This was repeated 10,000 times and the means for
F\textsubscript{ST} and D\textsubscript{XY} were taken across all
permutations. These permuted values were then compared to the empirical
mean F\textsubscript{ST} and D\textsubscript{XY} values for the
phenotype-associated SNPs using the \emph{t.test} function in R.

\subsection*{Analyses of recombination rate
variation:}

{\label{analyses-of-recombination-rate-variation}}

To test the effect of the recombination landscape on phenotype-genotype
association, we first estimated recombination rates using the published
Arctic charr linkage map (N=3,636) (Christensen et al., 2018) using
\emph{MareyMap v1.3.6} (Siberchicot et al., 2020). RAD loci from the
linkage map were aligned to the \emph{Salvelinus} sp. reference genome
with Bowtie2 (Langmead \& Salzberg, 2012) using the
\emph{-very-sensitive} setting. Loci were kept if they uniquely mapped
to one location, mapped to the same chromosome as all other loci on
their linkage group, and followed the orientation of the linkage map
(i.e., not reversed). The filtered dataset was used to estimate the
recombination rate across each chromosome using a spline algorithm. Spar
values were varied for each chromosome from 0.5 to 0.9, depending on
chromosome size, to best fit the data (Berloff et al., 2002).
Subsequently, \emph{WindowScanR v0.1} (available at:
\url{https://github.com/tavareshugo/WindowScanR}) was used to summarize
recombination rate values in 1 MB windows along the genome. All SNPs
were assigned to these windows using \emph{BEDtools}. A random subset of
100 SNPs was then selected and the mean recombination rate for those
SNPs was calculated based on their windows. This was repeated 10,000
times to generate a background mean recombination rate, which was then
compared to the mean recombination rate of the phenotype-associated SNPs
(based on their windows) using a t-test.

\section*{Results:}

{\label{results}}

\subsection*{Ecomorph divergence in head
shape:}

{\label{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 (\selectlanguage{greek}θ\selectlanguage{english}) was significantly different (P \textless{}
0.05) (Table 1). Similarly, almost all differences in trajectory lengths
between ecomorphs pairs (\selectlanguage{greek}Δ\selectlanguage{english}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 \textasciitilde{} ecomorph + lake + ecomorph x
lake) which found that the ecomorph x lake interaction term explained
most variation (\selectlanguage{greek}η\selectlanguage{english}\textsuperscript{2}\textsubscript{Eco*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).

\subsection*{Genomic regions associated with head
shape:}

{\label{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
(PVE\textsubscript{Head}) in the SNP dataset was 0.62 with the
proportion of phenotypic variance explained by large (``sparse'') effect
loci (PGE\textsubscript{Head}) 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 ([?]\textsubscript{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 \textless{} 0.05; Fig. 2A) with these SNPs broadly
distributed across the genome.

\subsection*{Genomic differentiation at SNPs associated with head
shape:}

{\label{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 (F\textsubscript{ST}) and absolute
divergence (D\textsubscript{XY}) between ecomorphs in lochs Dughaill and
Tay were elevated amongst associated SNPs when compared to the
background (Figure 4, 5) (Table S4). F\textsubscript{ST} and
D\textsubscript{XY} were significantly lower than background between
ecomorphs in Loch Awe for the associated SNPs. There was no significant
difference in F\textsubscript{ST} or D\textsubscript{XY} 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).

\subsection*{Genes and QTLs associated with head shape
variation:}

{\label{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 \selectlanguage{ngerman}±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 .

\subsection*{Ecomorph divergence in body
shape:}

{\label{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 (\selectlanguage{greek}Δ\selectlanguage{english}L) was significant with the
exception of the Awe-na Sealga comparison (Table 1). When comparing
angle of difference in trajectories (\selectlanguage{greek}θ\selectlanguage{english}), 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
(\selectlanguage{greek}η\selectlanguage{english}\textsuperscript{2}\textsubscript{Eco} = 0.301) for body shape
explained more variation than the ecomorph x lake term
(\selectlanguage{greek}η\selectlanguage{english}\textsuperscript{2}\textsubscript{Eco*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
(\selectlanguage{greek}η\selectlanguage{english}\textsuperscript{2}\textsubscript{Lake} = 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).

\subsection*{Genomic regions associated with body
shape:}

{\label{genomic-regions-associated-with-body-shape}}

A BSLMM found that the proportion of phenotypic variance in body shape
explained by the SNP dataset (PVE\textsubscript{Body}) was 0.82 and the
proportion of phenotypic variance explained by large (``sparse'') effect
loci (PGE\textsubscript{Body}) was 0.39 and the [?] (rho) value
([?]\textsubscript{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).

\subsection*{Genomic differentiation at SNPs associated with body
shape:}

{\label{genomic-differentiation-at-snps-associated-with-body-shape}}

For F\textsubscript{ST}, 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. D\textsubscript{XY} was significantly higher at Tay for the
associated SNPs while it was significantly lower at Loch Awe. There was
no significant difference in D\textsubscript{XY} 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).

\subsection*{Genes and QTLs associated with body
shape:}

{\label{genes-and-qtls-associated-with-body-shape}}

Relative to known genes in the \emph{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 .

\subsection*{Comparisons between head shape and body
shape:}

{\label{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
(\selectlanguage{greek}Δ\selectlanguage{english}L\textsubscript{Head} = 0.053 \selectlanguage{ngerman}± 0.035 s.d., mean
\selectlanguage{greek}θ\selectlanguage{english}\textsubscript{Head}= 69.87\selectlanguage{ngerman}° ± 17.54 s.d. mean) compared to body shape
(mean \selectlanguage{greek}Δ\selectlanguage{english}L\textsubscript{Body} = 0.016 \selectlanguage{ngerman}± 0.011 s.d., mean
\selectlanguage{greek}θ\selectlanguage{english}\textsubscript{Body}= 61.62\selectlanguage{ngerman}° ± 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\textsuperscript{-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
R\textsuperscript{2} = 0.99, P \textless{} 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 R\textsuperscript{2} = 0.29, P = 0.154) (Figure S6).

\section*{Discussion:}

{\label{discussion}}

Through our analyses, we identified genomic regions that underlie head
shape and/or body shape morphology in ecomorph pairs of Arctic charr. Of
these phenotype-associated SNPs, many (approximately one quarter: 50 of
212) were shared between head and body shape. The extent and direction
of divergence in head and body shape morphology were positively
correlated, suggesting a shared developmental basis for the two
phenotypes. The SNPs we found associated to each phenotype were often
found in genes related to morphology or anatomical development. Indeed,
independently previously identified QTLs in the genomic region of body
or head associated SNPs were often those related to morphology.

We found limited parallelism in shape morphology and in genomic
underpinnings with many population-specific patterns in the divergence
of head and body shape morphology between ecomorphs across pairs. Many
of the phenotype-associated SNPs were not present in all four pairs
likely due to polygenic genomic architectures and the incomplete
representation of the genome in our approach. The phenotype-associated
SNPs that were found to be shared were not highly diverged in all pairs
and did not appear to be under the same selective pressures. Body shape
in particular appears to be rather polygenic allowing for the genomic
underpinnings of the phenotypes to vary across lakes.

\subsection*{Limited parallelism in head and body shape divergences
across
replicates}

{\label{limited-parallelism-in-head-and-body-shape-divergences-across-replicates}}

While previous work on parallelism in these pairs suggested substantial
phenotypic parallelism in some linear traits related to head shape our
more sensitive geomorphometric approach to describe shape suggests
considerable phenotypic variation across replicated ecomorphs in
multivariate space for both head shape and body shape.

The lack of phenotypic parallelism between Loch Dughaill and the other
ecomorphs pairs may arise from its ecological distinctiveness compared
to the others. While benthivore morphs typically occupy the shallow
littoral zone of lakes as they do at the other three lakes in our study,
the benthivore morph at Loch Dughaill is a `profundal' benthivore and as
such utilises a much deeper part of the lake environment, and the lake
has a smaller area (Hooker et al., 2016; Jonsson \& Jonsson, 2001). The
divergence we see between the ecomorphs at Loch Dughaill, as indicated
by the PCA, is in line with other `profundal' charr morphs found across
the Holarctic, with the benthivore ecomorph having a deeper head and
body than the planktivore morph (Klemetsen, 2010; Skoglund et al.,
2015). This is the inverse pattern of the more common
benthivore-planktivore divergence where the benthivore morph has the
shallower, longer head and body shape that is seen in the other three
ecomorph pairs (Garduño-Paz et al., 2012). ~~

Evidence for the parallel evolution of head and body shape morphology
across ecomorphs has been shown in previous studies of Arctic charr (C.
E. Adams et al., 2008; Kristjánsson et al., 2012; Saltykova et al.,
2017). However, even disregarding the notably distinct Dughaill ecomorph
pair, the other three ecomorph pairs show limited evidence of
parallelism in shape morphology despite the parallel evolution of
ecomorphs themselves. Evolutionary divergence and thus phenotypic
trajectories are influenced by the interaction between environmental
variation and adaptive genetic variation. As a result, repeated ecomorph
divergences often have very different phenotypic trajectories for key
components of phenotypes, as seen in three-spine stickleback (Stuart et
al., 2017). Thus, the lack of strict parallelism seen in our study is
likely the result of known differences in the evolutionary histories of
these pairs and differing selective pressures in their local
environments, for example from differences in their ecosystems or diets
(Jacobs et al., 2020). ~Parallelism can also be generated by other
mechanisms, such as differences in gene expression, post-translational
modifications, and/or alternate splicing (Filteau et al., 2013; McGirr
\& Martin, 2018). Indeed, differences in splicing and gene expression
patterns showed parallel patterns across ecomorph pairs in a study on
three of the ecomorph pairs we investigated (Jacobs \& Elmer, 2021) and
provide explanations of alternative adaptive paths.

Of the two traits we tested, head shape had longer phenotypic trajectory
lengths and greater angles, suggesting that the ecomorphs are more
differentiated in head shape from one another and that head shape has
evolved in more distinct directions across lakes. Head shape is well
recognised as an important phenotype for foraging and prey
specialisation, while body shape is important for swimming behaviour and
habitat complexity in these and other fishes (C. E. Adams et al., 1998;
Skoglund et al., 2015; Webb, 1984). The considerable divergence in head
shape between ecomorphs within lakes suggests different prey
specialisations (Garduño-Paz et al. 2012; Hooker et al 2016) while their
body shapes and therefore perhaps swimming behaviours are more subtly
different. The high trajectory angles for head shape suggest notable
difference in foraging across lakes, whether that be due to the lake
environments or the actual species available as prey (Garduño-Paz \&
Adams, 2010). For both traits, Loch Tay and Dughaill showed notably
higher trajectory lengths than Awe and na Sealga and more evolutionary
divergence, also seen in the ecomorphs' genomic divergence
(F\textsubscript{ST} \textasciitilde{}1\% between ecomorphs in Awe and
na Sealga vs 9\% in Dughaill and Tay). ~This reflects likely what were
previously inferred to be recent sympatric divergences of ecomorph pairs
in Awe and na Sealga while Tay and Dughaill each have complex histories
of divergence and secondary contact between colonising lineages (Fenton
et al., 2023; Jacobs et al., 2020).

\subsection*{Genomic underpinnings of head and body shape across
lakes:}

{\label{genomic-underpinnings-of-head-and-body-shape-across-lakes}}

From our total of 212 SNPs that showed high associations with head
and/or body shape (Figure 2), we found more SNPs associated with body
shape than for head shape. Head shape was controlled by more
large-effect loci relative to body shape and may suggest that head shape
is controlled by fewer genes/pathways. In both cases these will be an
underestimate of actual associations because we have reduced
representation of the genome captured. We found for that for both head
and body shape only a small number of associated SNPs were diverged
between ecomorphs in all four pairs (Figure 3). This is line with what
has been suggested both in other Arctic charr studies and other salmonid
species, in which genetic differentiation between ecomorphs is largely
nonparallel across pairs bar at a few key genes (Jacobs et al. 2020;
Brachmann et al., 2022; Salisbury et al., 2020; Salisbury \& Ruzzante,
2021). Further to previous work, we found that the SNPs shared across
pairs were not highly differentiated between ecomorphs in all pairs
suggesting that while present, they are not critical to underlying the
phenotypic differences in each pair (Figure 4). These results also
suggest that the genomic underpinnings of each phenotype varies across
the lakes, likely contributing to the phenotypic differences we see
between pairs. The polygenic genomic underpinnings of both phenotypes,
as indicated by the numbers of associated SNPs identified, indicate that
there are multiple pathways that can achieve the same phenotypes hence
the lack of high divergence for the same SNPs across all lakes
(Salisbury \& Ruzzante, 2021).

Loch Dughaill and in particular Loch Tay often showed notable high
genomic divergence between ecomorphs for many of the associated SNPs for
each trait. Additionally, the associated SNPs for both traits showed
high D\textsubscript{XY} values compared to the background subsets at
both of these lakes. While increased levels of D\textsubscript{XY} or
F\textsubscript{ST} compared to genomic background can be indicative of
positive selection (Bamshad \& Wooding, 2003), they might also be
expected for loci resisting introgression following secondary contact
(Cruickshank \& Hahn, 2014), as is likely the case in Loch Tay and Loch
Dughaill (Jacobs et al., 2020). The associated SNPs we found are
widespread across the genome (Figure 2) indicating these are not single
linked regions of divergence as found in studies on Atlantic cod and
rainbow trout (Hemmer-Hansen et al., 2013; Pearse et al., 2019) but
instead are diffused and highly polygenic, similar to patterns for body
shape in lake whitefish (Laporte et al., 2015).

\subsection*{Functional genomic regions for head and body
shape:}

{\label{functional-genomic-regions-for-head-and-body-shape}}

Roughly half of the associated SNPs identified for each of head and body
shape were found within or proximal to an annotated gene in the charr
genome. A number of the GO terms that appeared as significantly
overrepresented or enriched in our study have been identified in other
studies investigating adaptive divergences or parallel evolution in
various fish species. Odontogenesis (GO:0042476), sensory perception of
sound (GO:0007605), blood vessel remodelling (GO:0001974), response to
muscle activity (GO:0014850), ventricular trabecula myocardium
morphogenesis (GO:0003222), common-partner SMAD protein phosphorylation
(GO:0007182), cellular response to ethanol (GO:0071361), and
neuromuscular synaptic transmission (GO:0007274) have shown significance
in other Arctic charr studies investigating ecomorph divergence
(Guðbrandsson et al., 2019; Salisbury et al., 2020). The GO terms for
associative learning (GO:0008306), regulation of cell shape
(GO:0008360), and UDP-glucuronate biosynthetic process (GO:0006065) also
appear in overrepresented groups in a study on the divergence of a
sympatric lake whitefish species pair (\emph{Corgeonus clupeaformis}) in
the USA (Hebert et al., 2013). Finally, in pupfishes (\emph{Cyprinodon.
sp.}), the divergent expression of a number of genes involved in cranial
skeletal system development was seen between different trophic
specialists with the GO term for this process (GO:1904888) significant
in our study (McGirr \& Martin, 2018). Differences in ossification rate
have been related to adaptive morphological differentiation in other
freshwater fish and the over enrichment or overexpression of genes
related to formation of various bones in our study indicates a similarly
important role in adaptive divergences between ecomorphs of Arctic charr
(Esin et al., 2018). Indeed, previous work has noted the importance of
differences in bone structure and sizes between different ecomorphs of
Arctic charr (Jónsdóttir et al., 2023; Kapralova et al., 2015). ~

The QTL database that we have developed allows us to explore whether
rapid replicated diversification of ecomorphs in different salmonid
species is underlined by the use of the same functional regions as has
been previously suggested for salinity tolerance (Jacobs et al., 2017;
Norman et al., 2012). Our results suggest that this is true to some
extent with QTLs related to body shape in lake trout and whitefish found
in proximity to SNPs that we identified as being associated with
phenotypic differences in Arctic charr ecomorphs. Whilst we only
identified a small number of QTLs located near the associated SNPs, this
is line with other work which suggests that shared basis for ecomorph
divergence across species may be limited (Salisbury \& Ruzzante, 2021).
This QTL marker database will be a valuable resource for future salmonid
research.

\section*{Conclusion:}

{\label{conclusion}}

Our results indicate differences in head and body shape responses to
ecological selection regimes across four replicate lakes. These
differing responses are likely enabled through the use of largely
different genetic bases across independent replicate ecomorph pairs.
Specifically, we found that only a small number of SNPs were shared
across all four pairs, suggesting limited genetic parallelism with these
shared SNPs under varying selective pressures across lakes. We found
that head and body shape morphology have a level of shared genetic
underpinnings in Arctic charr and that the genetics of these phenotypes
is shared to an extent across different salmonid species. Our analyses
highlight the complexity of the evolutionary genetics that underlie
parallel phenotypes across replicates. Further it we demonstrate the
power of using population replicates to resolve fundamental genetic and
evolutionary patterns from the noise of local variation.

\section*{Acknowledgements:}

{\label{acknowledgements}}

This research was supported by NatureScot (SF, CWB) and the School of
Biodiversity, One Health \& Veterinary Medicine at the University of
Glasgow (SF, CWB, CEA, KRE, AJ) and a Leverhulme Trust Early Career
Fellowship (ECF-2020-509) and a University of Glasgow Lord Kelvin/Adam
Smith Leadership Fellowship (AJ).

\section*{References:}

{\label{references}}

Adams, C. E., Fraser, D., Huntingford, F. A., Greer, R. B., Askew, C.
M., \& Walker, A. F. (1998). Trophic polymorphism among Arctic charr
from Loch Rannoch, Scotland. \emph{Journal of Fish Biology},
\emph{52}(6), 1259--1271. https://doi.org/10.1006/jfbi.1998.0676

Adams, C. E., \& Huntingford, F. A. (2002). Inherited differences in
head allometry in polymorphic Arctic charr from Loch Rannoch, Scotland.
\emph{Journal of Fish Biology}, \emph{60}(3), 515--520.
https://doi.org/10.1006/jfbi.2002.1867

Adams, C. E., Wilson, A. J., \& Ferguson, M. M. (2008). Parallel
divergence of sympatric genetic and body size forms of Arctic charr,
\emph{Salvelinus alpinus}, from two Scottish lakes. \emph{Biological
Journal of the Linnean Society}, \emph{95}(4), 748--757.
https://doi.org/10.1111/j.1095-8312.2008.01066.x

Adams, D. C., \& Collyer, M. L. (2013). Phenotypic trajectory analysis:
Comparison of shape change patterns in evolution and ecology.
\emph{Hystrix}, \emph{24}(1), 75--83.
https://doi.org/10.4404/hystrix-24.1-6298

Adams, D. C., \& Otárola-Castillo, E. (2013). Geomorph: An R package for
the collection and analysis of geometric morphometric shape data.
\emph{Methods in Ecology and Evolution}, \emph{4}(4), 393--399.
https://doi.org/10.1111/2041-210X.12035

Alexa, A., \& Rahnenfuhrer, J. (2020). \emph{topGO: Enrichment Analysis
for Gene Ontology} (v2.40.0).
https://bioc.ism.ac.jp/packages/3.13/bioc/manuals/topGO/man/topGO.pdf

Arendt, J., \& Reznick, D. (2008). Convergence and parallelism
reconsidered: what have we learned about the genetics of adaptation?
\emph{Trends in Ecology and Evolution}, \emph{23}(1), 26--32.
https://doi.org/10.1016/j.tree.2007.09.011

Baillie, S. M., Muir, A. M., Hansen, M. J., Krueger, C. C., \& Bentzen,
P. (2016). Genetic and phenotypic variation along an ecological gradient
in lake trout \emph{Salvelinus namaycush}. \emph{BMC Evolutionary
Biology}, \emph{16}(1). https://doi.org/10.1186/s12862-016-0788-8

Bamshad, M., \& Wooding, S. P. (2003). Signatures of natural selection
in the human genome. \emph{Nature Reviews Genetics}, \emph{4}(2),
99--111. https://doi.org/10.1038/nrg999

Beaujean, A. A. (2012). \emph{BaylorEdPsych: R Package for Baylor
University Educational Psychology Quantitative Courses.} (v0.5).
https://cran.r-project.org/src/contrib/Archive/BaylorEdPsych/

Berloff, N., Perola, M., \& Lange, K. (2002). Spline Methods for the
Comparison of Physical and Genetic Maps. \emph{Journal of Computational
Biology}, \emph{9}(3), 465--475.
https://doi.org/10.1089/106652702760138565

Bolnick, D. I., Barrett, R. D. H., Oke, K. B., Rennison, D. J., \&
Stuart, Y. E. (2018). (Non) Parallel Evolution. \emph{Annual Review of
Ecology, Evolution, and Systematics}, \emph{49}, 303--330.
https://doi.org/10.1146/annurev-ecolsys-110617-062240

Boulding, E. G., Culling, M., Glebe, B., Berg, P. R., Lien, S., \& Moen,
T. (2008). Conservation genomics of Atlantic salmon: SNPs associated
with QTLs for adaptive traits in parr from four trans-Atlantic
backcrosses. \emph{Heredity}, \emph{101}(4), 381--391.
https://doi.org/10.1038/hdy.2008.67

Brachmann, M. K., Parsons, K., Skúlason, S., Gaggiotti, O., \& Ferguson,
M. (2022). Variation in the genomic basis of parallel phenotypic and
ecological divergence in benthic and pelagic morphs of Icelandic Arctic
charr (\emph{Salvelinus alpinus}). \emph{Molecular Ecology},
\emph{31}(18), 4589--4599. https://doi.org/10.1111/mec.16625

Christensen, K. A., Rondeau, E. B., Minkley, D. R., Leong, J. S.,
Nugent, C. M., Danzmann, R. G., Ferguson, M. M., Stadnik, A., Devlin, R.
H., Muzzerall, R., Edwards, M., Davidson, W. S., \& Koop, B. F. (2018).
The Arctic charr (\emph{Salvelinus alpinus}) genome and transcriptome
assembly. \emph{PLoS ONE}, \emph{13}(9), 1--30.
https://doi.org/10.1371/journal.pone.0204076

Conte, G. L., Arnegard, M. E., Peichel, C. L., \& Schluter, D. (2012).
The probability of genetic parallelism and convergence in natural
populations. \emph{Proceedings of the Royal Society B: Biological
Sciences}, \emph{279}(1749), 5039--5047.
https://doi.org/10.1098/rspb.2012.2146

Cruickshank, T. E., \& Hahn, M. W. (2014). Reanalysis suggests that
genomic islands of speciation are due to reduced diversity, not reduced
gene flow. \emph{Molecular Ecology}, \emph{23}(13), 3133--3157.
https://doi.org/10.1111/mec.12796

Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E.,
DePristo, M. A., Handsaker, R. E., Lunter, G., Marth, G. T., Sherry, S.
T., McVean, G., \& Durbin, R. (2011). The variant call format and
VCFtools. \emph{Bioinformatics}, \emph{27}(15), 2156--2158.
https://doi.org/10.1093/bioinformatics/btr330

Doenz, C. J., Krähenbühl, A. K., Walker, J., Seehausen, O., \&
Brodersen, J. (2019). Ecological opportunity shapes a large Arctic charr
species radiation. \emph{Proceedings of the Royal Society B: Biological
Sciences}, \emph{286}(1913). https://doi.org/10.1098/rspb.2019.1992

Elmer, K. R. (2016). Genomic tools for new insights to variation,
adaptation, and evolution in the salmonid fishes: a perspective for
charr. \emph{Hydrobiologia}, \emph{783}(1), 191--208.
https://doi.org/10.1007/s10750-015-2614-5

Elmer, K. R., Fan, S., Kusche, H., Luise Spreitzer, M., Kautt, A. F.,
Franchini, P., \& Meyer, A. (2014). Parallel evolution of Nicaraguan
crater lake cichlid fishes via non-parallel routes. \emph{Nature
Communications}, \emph{5}. https://doi.org/10.1038/ncomms6168

Elmer, K. R., \& Meyer, A. (2011). Adaptation in the age of ecological
genomics: Insights from parallelism and convergence. \emph{Trends in
Ecology \& Evolution}, \emph{26}, 298--306.
https://doi.org/10.1016/j.tree.2011.02.008

Esin, E. v., Markevich, G. N., \& Pichugin, M. Y. (2018). Juvenile
divergence in adaptive traits among seven sympatric fish ecomorphs
arises before moving to different lacustrine habitats. \emph{Journal of
Evolutionary Biology}, \emph{31}(7), 1018--1034.
https://doi.org/10.1111/jeb.13283

Fang, B., Kemppainen, P., Momigliano, P., Feng, X., \& Merilä, J.
(2020). On the causes of geographically heterogeneous parallel evolution
in sticklebacks. \emph{Nature Ecology and Evolution}, \emph{4},
1105--1115. https://doi.org/10.1038/s41559-020-1222-6

Fenton, S., Elmer, K. R., Bean, C. W., \& Adams, C. E. (2023). How
glaciation impacted evolutionary history and contemporary genetic
diversity of flora and fauna in the British Isles. \emph{Scottish
Geographical Journal}, 1--21.
https://doi.org/10.1080/14702541.2023.2231407

Filteau, M., Pavey, S. A., St-Cyr, J., \& Bernatchez, L. (2013). Gene
coexpression networks reveal key drivers of phenotypic divergence in
lake whitefish. \emph{Molecular Biology and Evolution}, \emph{30}(6),
1384--1396. https://doi.org/10.1093/molbev/mst053

Garduño-Paz, M. V., \& Adams, C. E. (2010). Discrete prey availability
promotes foraging segregation and early divergence in Arctic charr,
Salvelinus alpinus. \emph{Hydrobiologia}, \emph{650}(1), 15--26.
https://doi.org/10.1007/s10750-009-0055-8

Garduño-Paz, M. V., Adams, C. E., Verspoor, E., Knox, D., \& Harrod, C.
(2012). Convergent evolutionary processes driven by foraging opportunity
in two sympatric morph pairs of Arctic charr with contrasting
post-glacial origins. \emph{Biological Journal of the Linnean Society},
\emph{106}(4), 794--806.
https://doi.org/10.1111/j.1095-8312.2012.01906.x

Guðbrandsson, J., Kapralova, K. H., Franzdóttir, S. R.,
Bergsveinsdóttir, Þ. M., Hafstað, V., Jónsson, Z. O., Snorrason, S. S.,
\& Pálsson, A. (2019). Extensive genetic differentiation between
recently evolved sympatric Arctic charr morphs. \emph{Ecology and
Evolution}, \emph{9}(19), 10964--10983.
https://doi.org/10.1002/ece3.5516

Hebert, F. O., Renaut, S., \& Bernatchez, L. (2013). Targeted sequence
capture and resequencing implies a predominant role of regulatory
regions in the divergence of a sympatric lake whitefish species pair
(\emph{Coregonus clupeaformis}). \emph{Molecular Ecology},
\emph{22}(19), 4896--4914. https://doi.org/10.1111/mec.12447

Hemmer-Hansen, J., Nielsen, E. E., Therkildsen, N. O., Taylor, M. I.,
Ogden, R., Geffen, A. J., Bekkevold, D., Helyar, S., Pampoulie, C.,
Johansen, T., \& Carvalho, G. R. (2013). A genomic island linked to
ecotype divergence in Atlantic cod. \emph{Molecular Ecology},
\emph{22}(10), 2653--2667. https://doi.org/10.1111/mec.12284

Hooker, O. E., Barry, J., Van Leeuwen, T. E., Lyle, A., Newton, J.,
Cunningham, P., \& Adams, C. E. (2016). Morphological, ecological and
behavioural differentiation of sympatric profundal and pelagic Arctic
charr (\emph{Salvelinus alpinus}) in Loch Dughaill Scotland.
\emph{Hydrobiologia}, \emph{783}(1), 209--221.
https://doi.org/10.1007/s10750-015-2599-0

Jacobs, A., Carruthers, M., Yurchenko, A., Gordeeva, N. V., Alekseyev,
S. S., Hooker, O., Leong, J. S., Minkley, D. R., Rondeau, E. B., Koop,
B. F., Adams, C. E., \& Elmer, K. R. (2020). Parallelism in
eco-morphology and gene expression despite variable evolutionary and
genomic backgrounds in a Holarctic fish. \emph{PLOS Genetics},
\emph{16}(4), e1008658. https://doi.org/10.1371/journal.pgen.1008658

Jacobs, A., \& Elmer, K. R. (2021). Alternative splicing and gene
expression play contrasting roles in the parallel phenotypic evolution
of a salmonid fish. \emph{Molecular Ecology}, \emph{November 2020},
1--15. https://doi.org/10.1111/mec.15817

Jacobs, A., Womack, R., Chen, M., Gharbi, K., \& Elmer, K. R. (2017).
Significant synteny and colocalization of ecologically relevant
quantitative trait loci within and across species of salmonid fishes.
\emph{Genetics}, \emph{207}(2), 741--754.
https://doi.org/10.1534/genetics.117.300093

Jensen, H., Kiljunen, M., Knudsen, R., \& Amundsen, P. A. (2017).
Resource partitioning in food, space and time between Arctic charr
(\emph{Salvelinus alpinus}), brown trout (\emph{Salmo trutta}) and
european whitefish (\emph{Coregonus lavaretus}) at the southern edge of
their continuous coexistence. \emph{PLoS ONE}, \emph{12}(1), 1--18.
https://doi.org/10.1371/journal.pone.0170582

Jónsdóttir, G. Ó., von Elm, L.-M., Ingimarsson, F., Tersigni, S.,
Snoarrason, S. S., Pálsson, A., \& Steele, S. E. (2023). Diversity in
the internal functional feeding elements of sympatric morphs of Arctic
charr (\emph{Salvelinus alpinus}). \emph{BioRxiv}.
https://doi.org/10.1101/2023.02.17.528955

Jonsson, B., \& Jonsson, N. (2001). Polymorphism and speciation in
Arctic charr. \emph{Journal of Fish Biology}, \emph{58}(3), 605--638.
https://doi.org/10.1006/jfbi.2000.1515

Kapralova, K. H., Ias, Z., Onsson, O. J., Palsson, A., Idur, S., Le
Deuff, S., Kristjánsson, B. K., \& Snorrason, S. S. (2015). Bones in
Motion: Ontogeny of Craniofacial Development in Sympatric Arctic Charr
Morphs. \emph{Developmental Dynamics}, \emph{244}, 1168--1178.
https://doi.org/10.1002/dvdy

Klemetsen, A. (2010). The Charr Problem Revisited: Exceptional
Phenotypic Plasticity Promotes Ecological Speciation in Postglacial
Lakes. \emph{Freshwater Reviews}, \emph{3}(1), 49--74.
https://doi.org/10.4290/frj-3.1.3

Kowalko, J. E., Rohner, N., Linden, T. A., Rompani, S. B., Warren, W.
C., Borowsky, R., Tabin, C. J., Jeffery, W. R., \& Yoshizawa, M. (2013).
Convergence in feeding posture occurs through different genetic loci in
independently evolved cave populations of \emph{Astyanax mexicanus}.
\emph{Proceedings of the National Academy of Sciences}, \emph{110}(42),
16933--16938. https://doi.org/10.1073/pnas.1317192110

Kristjánsson, B. K., Skúlason, S., Snorrason, S. S., \& Noakes, D. L. G.
(2012). Fine-scale parallel patterns in diversity of small benthic
Arctic charr (\emph{Salvelinus alpinus}) in relation to the ecology of
lava/groundwater habitats. \emph{Ecology and Evolution}, \emph{2}(6),
1099--1112. https://doi.org/10.1002/ece3.235

Küttner, E., Parsons, K. J., Easton, A. A., Skúlason, S., Danzmann, R.
G., \& Ferguson, M. M. (2014). Hidden genetic variation evolves with
ecological specialization: The genetic basis of phenotypic plasticity in
Arctic charr ecomorphs. \emph{Evolution and Development}, \emph{16}(4),
247--257. https://doi.org/10.1111/ede.12087

Landry, L., Vincent, W. F., \& Bernatchez, L. (2007). Parallel evolution
of lake whitefish dwarf ecotypes in association with limnological
features of their adaptive landscape. \emph{Journal of Evolutionary
Biology}, \emph{20}(3), 971--984.
https://doi.org/10.1111/j.1420-9101.2007.01304.x

Langerhans, R. B., \& DeWitt, T. J. (2004). Shared and unique features
of evolutionary diversification. \emph{American Naturalist},
\emph{164}(3), 335--349. https://doi.org/10.1086/422857

Langmead, B., \& Salzberg, S. L. (2012). Fast gapped-read alignment with
Bowtie 2. \emph{Nature Methods}, \emph{9}(4), 357--359.
https://doi.org/10.1038/nmeth.1923

Laporte, M., Rogers, S. M., Dion-Côté, A. M., Normandeau, E., Gagnaire,
P. A., Dalziel, A. C., Chebib, J., \& Bernatchez, L. (2015). RAD-QTL
mapping reveals both genome-level parallelism and different genetic
architecture underlying the evolution of body shape in lake Whitefish
(\emph{Coregonus clupeaformis}) species pairs. \emph{G3: Genes, Genomes,
Genetics}, \emph{5}(7), 1481--1491.
https://doi.org/10.1534/g3.115.019067

Larson, W. A., Dann, T. H., Limborg, M. T., McKinney, G. J., Seeb, J.
E., \& Seeb, L. W. (2019). Parallel signatures of selection at genomic
islands of divergence and the major histocompatibility complex in
ecotypes of sockeye salmon across Alaska. \emph{Molecular Ecology},
\emph{28}(9), 2254--2271. https://doi.org/10.1111/mec.15082

Magalhaes, I. S., Whiting, J. R., D'Agostino, D., Hohenlohe, P. A.,
Mahmud, M., Bell, M. A., Skúlason, S., \& MacColl, A. D. C. (2021).
Intercontinental genomic parallelism in multiple three-spined
stickleback adaptive radiations. \emph{Nature Ecology and Evolution},
\emph{5}(2), 251--261. https://doi.org/10.1038/s41559-020-01341-8

McGirr, J. A., \& Martin, C. H. (2018). Parallel evolution of gene
expression between trophic specialists despite divergent genotypes and
morphologies. \emph{Evolution Letters}, \emph{2}(2), 62--75.
https://doi.org/10.1002/evl3.41

Money, D., Gardner, K., Migicovsky, Z., Schwaninger, H., Zhong, G. Y.,
\& Myles, S. (2015). LinkImpute: Fast and accurate genotype imputation
for nonmodel organisms. \emph{G3: Genes, Genomes, Genetics},
\emph{5}(11), 2383--2390. https://doi.org/10.1534/g3.115.021667

Norman, J. D., Danzmann, R. G., Glebe, B., \& Ferguson, M. M. (2011).
The genetic basis of salinity tolerance traits in Arctic charr
(\emph{Salvelinus alpinus}). \emph{BMC Genetics}, \emph{12}, 1--12.
https://doi.org/10.1186/1471-2156-12-81

Norman, J. D., Robinson, M., Glebe, B., Ferguson, M. M., \& Danzmann, R.
G. (2012). Genomic arrangement of salinity tolerance QTLs in salmonids:
A comparative analysis of Atlantic salmon (\emph{Salmo salar}) with
Arctic charr (\emph{Salvelinus alpinus}) and rainbow trout
(\emph{Oncorhynchus mykiss}). \emph{BMC Genomics}, \emph{13}(1).
https://doi.org/10.1186/1471-2164-13-420

Oke, K. B., Rolshausen, G., LeBlond, C., \& Hendry, A. P. (2017). How
parallel is parallel evolution? A comparative analysis in fishes.
\emph{American Naturalist}, \emph{190}(1).
https://doi.org/10.1086/691989

Pearse, D. E., Barson, N. J., Nome, T., Gao, G., Campbell, M. A.,
Abadía-Cardoso, A., Anderson, E. C., Rundio, D. E., Williams, T. H.,
Naish, K. A., Moen, T., Liu, S., Kent, M., Moser, M., Minkley, D. R.,
Rondeau, E. B., Brieuc, M. S. O., Sandve, S. R., Miller, M. R., \ldots{}
Lien, S. (2019). Sex-dependent dominance maintains migration supergene
in rainbow trout. \emph{Nature Ecology and Evolution}, \emph{3}(12),
1731--1742. https://doi.org/10.1038/s41559-019-1044-6

Quinlan, A. R., \& Hall, I. M. (2010). BEDTools: A flexible suite of
utilities for comparing genomic features. \emph{Bioinformatics},
\emph{26}(6), 841--842. https://doi.org/10.1093/bioinformatics/btq033

R Core Team. (2021). \emph{R: A language and environment for statistical
computing}. R Foundation for Statistical Computing.
https://www.R-project.org/

Ravinet, M., Faria, R., Butlin, R. K., Galindo, J., Bierne, N.,
Rafajlovi\selectlanguage{polish}ć, \selectlanguage{english}M., Noor, M. A. F., Mehlig, B., \& Westram, A. M. (2017).
Interpreting the genomic landscape of speciation: a road map for finding
barriers to gene flow. \emph{Journal of Evolutionary Biology},
\emph{30}(8), 1450--1477. https://doi.org/10.1111/jeb.13047

Ravinet, M., Westram, A., Johannesson, K., Butlin, R., Andr\selectlanguage{ngerman}é, C., \&
Panova, M. (2016). Shared and nonshared genomic divergence in parallel
ecotypes of \emph{Littorina saxatilis} at a local scale. \emph{Molecular
Ecology}, \emph{25}(1), 287--305. https://doi.org/10.1111/mec.13332

Rochette, N. C., Rivera-Colón, A. G., \& Catchen, J. M. (2019). Stacks
2: Analytical methods for paired-end sequencing improve RADseq-based
population genomics. \emph{Molecular Ecology}, \emph{28}(21),
4737--4754. https://doi.org/10.1111/mec.15253

Salisbury, S. J., McCracken, G. R., Perry, R., Keefe, D., Layton, K. K.
S., Kess, T., Nugent, C. M., Leong, J. S., Bradbury, I. R., Koop, B. F.,
Ferguson, M. M., \& Ruzzante, D. E. (2020). Limited genetic parallelism
underlies recent, repeated incipient speciation in geographically
proximate populations of an Arctic fish (\emph{Salvelinus alpinus}).
\emph{Molecular Ecology}, \emph{29}(22), 4280--4294.
https://doi.org/10.1111/mec.15634

Salisbury, S. J., \& Ruzzante, D. E. (2021). Genetics causes and
Consequences of Sympatric Morph Divergence in \emph{Salmonidae}: A
Search for Mechanisms. \emph{Annual Review of Animal Biosciences}.
https://doi.org/10.1146/annurev-animal-021419-

Saltykova, E., Siwertsson, A., \& Knudsen, R. (2017). Parallel
phenotypic evolution of skull-bone structures and head measurements of
Arctic charr morphs in two subarctic lakes. \emph{Environmental Biology
of Fishes}, \emph{100}(2), 137--148.
https://doi.org/10.1007/s10641-016-0564-z

Schluter, D. (1996). Ecological speciation in postglacial fishes.
\emph{Philosophical Transactions of the Royal Society B: Biological
Sciences}, \emph{351}(1341), 807--814.
https://doi.org/10.1098/rstb.1996.0075

Schneider, C. A., Rasband, W. S., \& Eliceiri, K. W. (2012). NIH Image
to ImageJ: 25 years of image analysis. \emph{Nature Methods},
\emph{9}(7), 671--675. https://doi.org/10.1038/nmeth.2089

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage,
D., Amin, N., Schwikowski, B., \& Ideker, T. (2003). Cytoscape: A
software Environment for integrated models of biomolecular interaction
networks. \emph{Genome Research}, \emph{13}(11), 2498--2504.
https://doi.org/10.1101/gr.1239303

Siberchicot, A., Rezvoy, C., Charif, D., Gueguen, L., \& Marais, G.
(2020). \emph{MareyMap:Estimation of Meiotic Recombination Rates Using
Marey Maps} (v1.3.6). https://github.com/aursiber/MareyMap

Siwertsson, A., Knudsen, R., Adams, C. E., Præbel, K., \& Amundsen, P.
A. (2013). Parallel and non-parallel morphological divergence among
foraging specialists in European whitefish (\emph{Coregonus lavaretus}).
\emph{Ecology and Evolution}, \emph{3}(6), 1590--1602.
https://doi.org/10.1002/ece3.562

Skoglund, S., Siwertsson, A., Amundsen, P. A., \& Knudsen, R. (2015).
Morphological divergence between three Arctic charr morphs - the
significance of the deep-water environment. \emph{Ecology and
Evolution}, \emph{5}(15), 3114--3129. https://doi.org/10.1002/ece3.1573

Smith, S. R., Amish, S. J., Bernatchez, L., le Luyer, J., Wilson, C.,
Boeberitz, O., Luikart, G., \& Scribner, K. T. (2020). Mapping of
Adaptive Traits Enabled by a High-Density Linkage Map for Lake Trout.
\emph{G3: Genes, Genomes, Genetics}, \emph{10}(June), g3.401184.2020.
https://doi.org/10.1534/g3.120.401184

Stuart, Y. E., Veen, T., Weber, J. N., Hanson, D., Ravinet, M., Lohman,
B. K., Thompson, C. J., Tasneem, T., Doggett, A., Izen, R., Ahmed, N.,
Barrett, R. D. H., Hendry, A. P., Peichel, C. L., \& Bolnick, D. I.
(2017). Contrasting effects of environment and genetics generate a
continuum of parallel evolution. \emph{Nature Ecology and Evolution},
\emph{1}(6), 1--7. https://doi.org/10.1038/s41559-017-0158

Supek, F., Bo\selectlanguage{polish}š\selectlanguage{english}njak, M., \selectlanguage{polish}Š\selectlanguage{english}kunca, N., \& \selectlanguage{polish}Š\selectlanguage{english}muc, T. (2011). Revigo
summarizes and visualizes long lists of gene ontology terms. \emph{PLoS
ONE}, \emph{6}(7). https://doi.org/10.1371/journal.pone.0021800

Webb, P. W. (1984). Body Form, Locomotion and Foraging in Aquatic
Vertebrates. \emph{American Zoologist}, \emph{24}, 107--120.
https://academic.oup.com/icb/article/24/1/107/104211

Weinstein, S. Y., Thrower, F. P., Nichols, K. M., \& Hale, M. C. (2019).
A large-scale chromosomal inversion is not associated with life history
development in rainbow trout from Southeast Alaska. \emph{PLoS ONE},
\emph{14}(9). https://doi.org/10.1371/journal.pone.0223018

Wickham, H. (2016). \emph{ggplot2: Elegant Graphics for Data Analysis}.
Springer-Verlag New York. https://ggplot2.tidyverse.org

Zheng, X., Levine, D., Shen, J., Gogarten, S. M., Laurie, C., \& Weir,
B. S. (2012). A high-performance computing toolset for relatedness and
principal component analysis of SNP data. \emph{Bioinformatics},
\emph{28}(24), 3326--3328. https://doi.org/10.1093/bioinformatics/bts606

Zhou, X., Carbonetto, P., \& Stephens, M. (2013). Polygenic Modeling
with Bayesian Sparse Linear Mixed Models. \emph{PLoS Genetics},
\emph{9}(2). https://doi.org/10.1371/journal.pgen.1003264

Zhou, X., \& Stephens, M. (2012). Genome-wide efficient mixed-model
analysis for association studies. \emph{Nature Genetics}, \emph{44}(7),
821--824. https://doi.org/10.1038/ng.2310

Zhou, X., \& Stephens, M. (2014). Efficient multivariate linear mixed
model algorithms for genome-wide association studies. \emph{Nature
Methods}, \emph{11}(4), 407--409. https://doi.org/10.1038/nmeth.2848

\section*{Data Accessibility and
Benefit-Sharing:}

{\label{data-accessibility-and-benefit-sharing}}

The SNP VCF file, landmark data, and the QTL database are available on
the University of Glasgow's data repository, Enlighten (DOI:
http://dx.doi.org/10.5525/gla.researchdata.1502 {[}closed until
acceptance{]}). Benefits from this research accrue from the sharing of
our data and results on public databases.

\section*{Author Contributions:}

{\label{author-contributions}}

KRE and AJ conceived the study with input from CB and CA. Read mapping
and landmarking was done by AJ. Data analyses were performed by SF with
analytics support from AJ. SF wrote the original draft with support from
KRE and all authors contributed to subsequent versions of the
manuscript.

\section*{Disclosure statement:}

{\label{disclosure-statement}}

The authors report there are no competing interests to declare.

\subsection*{Table 1:}

{\label{table-1}}

Table 1. Phenotypic trajectory analysis comparisons across ecomorph
pairs for head shape and body shape. The difference in trajectory
length, the magnitude of change, between ecomorph pair is indicated by
\selectlanguage{greek}Δ\selectlanguage{english}L. The angle between trajectories is indicated as \selectlanguage{greek}θ. \selectlanguage{english}The significance
values are provided for each comparison.\selectlanguage{english}
\begin{longtable}[]{@{}llllll@{}}
\toprule
\textbf{Phenotype} & \textbf{Ecomorph pair comparison} & \selectlanguage{greek}\textbf{ΔL
(magnitude of change)}\selectlanguage{english} & \textbf{p value} & \selectlanguage{greek}\textbf{θ (direction of
change)}\selectlanguage{english} & \textbf{p value}\tabularnewline
\midrule
\endhead
\textbf{Head shape} & Awe-Dughaill & 0.0347 & 0.003 & 67.51\selectlanguage{ngerman}° &
0.001\tabularnewline
& Awe-naSealga & 0.0032 & 0.777 & 65.74\selectlanguage{ngerman}° & 0.002\tabularnewline
& Awe-Tay & 0.0901 & 0.001 & 57.29\selectlanguage{ngerman}° & 0.007\tabularnewline
& Dughaill-naSealga & 0.0379 & 0.002 & 78.37\selectlanguage{ngerman}° & 0.001\tabularnewline
& Dughaill-Tay & 0.0554 & 0.001 & 100.56\selectlanguage{ngerman}° & 0.001\tabularnewline
& naSealga-Tay & 0.0933 & 0.001 & 47.88\selectlanguage{ngerman}° & 0.027\tabularnewline
\textbf{Body shape} & Awe-Dughaill & 0.0097 & 0.004 & 92.48\selectlanguage{ngerman}° &
0.001\tabularnewline
& Awe-naSealga & 0.0013 & 0.726 & 55.27\selectlanguage{ngerman}° & 0.004\tabularnewline
& Awe-Tay & 0.028 & 0.001 & 39.64\selectlanguage{ngerman}° & 0.038\tabularnewline
& Dughaill-naSealga & 0.011 & 0.001 & 70.65\selectlanguage{ngerman}° & 0.001\tabularnewline
& Dughaill-Tay & 0.0183 & 0.001 & 79.24\selectlanguage{ngerman}° & 0.001\tabularnewline
& naSealga-Tay & 0.0293 & 0.001 & 34.20\selectlanguage{ngerman}° & 0.118\tabularnewline
\bottomrule
\end{longtable}

\subsection*{Table 2:}

{\label{table-2}}

Table 2: Table of associated SNPs found within \selectlanguage{ngerman}±100kbp of salmonid QTLs
mapped to the \emph{Salvelinus sp.} genome. Which phenotype the SNP is
associated with, its position, and position of QTL marker in question
are all indicated. The QTL type and species of origin are indicated. QTL
name refers to the designated the QTL was given in the whole QTL
database found in Table S1.\selectlanguage{english}
\begin{longtable}[]{@{}llllllll@{}}
\toprule
\textbf{SNP phenotype} & \textbf{SNP phenotype} & \textbf{QTL species} &
\textbf{QTL marker} & \textbf{QTL type} & \textbf{Chromosome} &
\textbf{SNP position} & \textbf{QTL position}\tabularnewline
\midrule
\endhead
Head and Body & \emph{C.clupeaformis} & \emph{C.clupeaformis} &
Cocl\_BS\_096 & Body shape & NC\_036838.1 & 46764741 &
46685809\tabularnewline
Head and Body & \emph{S.namaycush} & \emph{S.namaycush} & Sna\_BS\_069 &
Body shape & NW\_019942687.1 & 358354 & 361478\tabularnewline
Head and Body & \emph{S.namaycush} & \emph{S.namaycush} & Sna\_BS\_069 &
Body shape & NW\_019942687.1 & 358355 & 361478\tabularnewline
Body & \emph{S.alpinus} & \emph{S.alpinus} & Sal\_BW\_053 & Body weight
& NC\_036871.1 & 32030838 & 31993282\tabularnewline
Body & \emph{S.namaycush} & \emph{S.namaycush} & Sna\_BS\_092 & Body
shape & NC\_036854.1 & 21330886 & 21293280\tabularnewline
\bottomrule
\end{longtable}

\subsection*{Figure 1:}

{\label{figure-1}}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

Figure 1. Sample locations and morphological analysis. (A) Map of
Scotland showing the sampling locations. (B) Pictures showing an example
ecomorph pair, here from Loch Tay. Principal component analysis for
landmark analysis for Head shape (C) and Body shape (D). Individuals are
coloured by lake with larger points representing the mean values for
each ecomorph with a line connecting means in each pair. Vector plots
display how head shape (E) and body shape (F) morphology changes across
PC1 with diagrams showing where each landmark was placed on the fish for
each phenotype. Grey points show landmark positions at the minimum PC1
score and vectors show how landmark positions changes at the maximum PC1
score.

\subsection*{Figure 2:}

{\label{figure-2}}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

Figure 2: Manhattan plots for genomic location of SNPs highly associated
with morphology. A) Shown are 67 SNPs associated with head shape across
the four lakes of two ecomorphs of charr. B) Shown are 144 SNPs
associated with body shape. SNPs associated with head or body shape are
highlighted in blue; SNPs associated with both head shape and body shape
are highlighted in orange (38 SNPs). Red asterisk indicate SNPs shared
across all four ecomorph pairs for head shape (N=2) and body shape
(N=9). SNPs on unanchored scaffolds are not pictured.

\subsection*{Figure 3}

{\label{figure-3}}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

Figure 3: Venn diagram of SNPs associated with head shape (A) and body
shape (B) and how they are shared across each combination of different
lake pairs.

\subsection*{Figure 4:}

{\label{figure-4}}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image7/image7}

\end{center}
\end{figure}

Figure 4: Genetic differentiation (F\textsubscript{ST}) between
ecomorphs at each lake for different associated SNPs datasets. Body
refers to all SNPs associated for body shape and Head refers to all SNPs
associated with head shape. Body-shared dataset is just the body shape
SNPs found in all four ecomorph pairs and Head-shared is the equivalent
for head shape SNPs. Background SNPs refers to a randomly selected
background subset of SNPs used for comparisons. Red diamonds are the
mean value for that dataset. Significance of difference in means is
indicated by NS (P \textgreater{} 0.05), * (P \textless{} 0.05), ** (P
\textless{} 0.01), *** (P \textless{} 0.001).

\subsection*{Figure 5:}

{\label{figure-5}}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image9/image9}

\end{center}
\end{figure}

Figure 5: Absolute divergence (D\textsubscript{XY)} between ecomorphs at
each lake for different associated SNPs datasets. Body refers to all
SNPs associated for body shape and Head refers to all SNPs associated
with head shape. Body-shared dataset is just the body shape SNPs found
in all four ecomorph pairs and Head-shared is the equivalent for head
shape SNPs. Background SNPs refers to a randomly selected background
subset of SNPs used for comparisons. Red diamonds are the mean value for
that dataset. Significance of difference in means is indicated by NS (P
\textgreater{} 0.05), * (P \textless{} 0.05), ** (P \textless{} 0.01),
*** (P \textless{} 0.001).

\subsection*{Figure 6:}

{\label{figure-6}}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image11/image11}

\end{center}
\end{figure}

Figure 6. Comparing differences in the magnitude of phenotypic change
between ecomorphs for head (\selectlanguage{greek}Δ\selectlanguage{english}Lhead) and body shape
(\selectlanguage{greek}Δ\selectlanguage{english}L\textsubscript{Body}) across ecomorph pairs. Each ecomorph pair
comparison is indicated on the plot. Dug refers to Dughaill and Sea to
na Sealga.

\selectlanguage{english}
\FloatBarrier
\end{document}

