Further evidence from common garden rearing experiments of heritable traits separating lean and siscowet lake charr (Salvelinus namaycush) ecotypes.

Genetic evidence of selection for complex and polygenically regulated phenotypes can easily become masked by neutral population genetic structure and phenotypic plasticity. Without direct evidence of genotype-phenotype associations it can be difficult to conclude to what degree a phenotype is heritable or a product of environment. Common garden laboratory studies control for environmental stochasticity and help to determine the mechanism that regulate traits. Here we assess lipid content, growth, weight, and length variation in full and hybrid F1 crosses of deep and shallow water sympatric lake charr ecotypes reared for nine years in a common garden experiment. Redundancy analysis (RDA) and quantitative-trait-loci (QTL) genomic scans are used to identify associations between genotypes at 19,714 single nucleotide polymorphisms (SNPs) aligned to the lake charr genome and individual phenotypes to determine the role that genetic inheritance plays in ecotype phenotypic diversity. Lipid content, growth, length, and weight differed significantly among lake charr crosses throughout the experiment suggesting that pedigree plays a large roll in lake charr development. Polygenic scores of 15 SNPs putatively associated with lipid content and/or condition factor indicated that ecotype distinguishing traits are polygenically regulated and additive. A QTL identified on chromosome 38 contained >200 genes, some of which were associated with lipid metabolism and growth, demonstrating the complex nature of ecotype diversity. The results of our common garden study further indicate that lake charr ecotypes observed in nature are predetermined at birth and that ecotypes differ fundamentally in lipid metabolism and growth.

INTRODUCTION

Strong ecological gradients and intraspecific competition can cause species to form sympatric morphs or ecotypes (Lowry, 2012; Martin & Phennig, 2010). Ecological gradients can create habitat heterogeneity and alternate niches that allow species to establish resource polymorphisms that avoid competition (MacArthur, 1958; Skúlason & Smith, 1995). Understanding how intraspecific morphological divergence evolves along gradients, and to what degree phenotypic traits are environmentally and genetically regulated, continues to be an important question in evolutionary biology (Martin & Phennig, 2010; Skúlason et al., 2019). Lack of clarity about the relative role that environment and genetics have in determining phenotype can influence conclusions about the origin of intraspecific variation and the ecological forces driving selection (Skúlason et al., 2019). If phenotypic traits are genetically constrained, their maintenance may depend on there being limited gene flow between sympatric populations. At the same time, environmentally regulated phenotypes may persist even with gene flow. Therefore, determining the relative contributions of environment and genetics on phenotype is an important aspect of understanding the evolution of sympatric diversity.

Depth is a strong ecological gradient that commonly influences aquatic and marine environments. Light, temperature, and pressure change rapidly with depth, changing resource availability, metabolic processes, and foraging strategies (Farré et al., 2016). Because of this, aquatic species are often adapted to a specific depth or niche within the water column (Friedman et al., 2020). Deep‐water benthic species may have eyes adapted to low light levels, a different lipid composition or amount for maintaining buoyancy, and dark pigmentation to camouflage from visual predators coming from above (Farré et al., 2016; Radnaeva et al., 2017; Phleger, 1998; Yerlikaya et al., 2013). In contrast, pelagic fish that spend more time in open waters closer to the surface tend to be more fusiform, have smaller eyes, lighter coloration and use a swim bladder for buoyancy (Friedman et al., 2020; Pelster, 2021). Thus, depth is an ecological gradient that has long been recognized as a strong selective force contributing to the high diversity of fishes.

Similar patterns of depth and habitat‐related fish diversity have been observed in freshwater lakes (Eadie & Keast, 1984). In some cases, fish species in northern lakes express resource polymorphisms and fill multiple niches within their environment (e.g. Parsons & Robinson, 2007; Parsons et al., 2010). These polymorphisms are most common in the Salmonidae (trout, salmon, whitefish, charr), Gasterosteidae (sticklebacks) and Osmeridae (smelts) but also include fish in other families (Taylor, 1999). Under stable conditions, divergent selection helps to maintain sympatric phenotype differences and eventually lead to the formation of genetically distinct ecotypes (Linck et al., 2020). However, if the fitness benefits that maintain distinction weaken, renewed gene flow may begin to homogenize a population (Johansson & Ripa, 2006; Taylor et al., 2006; Thompson et al., 2019). Especially in lakes, sympatric ecotypes often lack physical reproductive barriers, and therefore reproductive isolation may be dependent on divergent selection.

Phenotypic differences between sympatric ecotypes can appear to be almost entirely plastic within each generation (Robinson & Wilson, 1996) or almost entirely heritable (Nosil et al., 2003). However, in many cases the mechanism behind phenotypic diversity is much less clear. Lake charr (Salvelinus namaycush) exhibit a great deal of phenotypic and ecological variation between and within lacustrine systems in North America (Chavarie et al., 2021). Many lakes contain sympatric lake charr ecotypes that can be distinguished by morphology, habitat, and diet (e.g. Baillie et al., 2018; Moore & Bronte, 2001). Most genetic studies indicate that lake charr ecotypes are not reproductively isolated, suggesting that phenotype may be environmentally controlled (Bernatchez et al., 2016; Chavarie et al., 2017). However, others have found evidence that differences among lake charr ecotypes is genetically regulated (Kissinger et al., 2018; Perreault‐Payette et al., 2017), even when gene flow among sympatric ecotypes is still prevalent.

The Laurentian Great Lakes of North America once held numerous distinct ecotypes of lake charr; however, overharvest and invasive lamprey predation are believed to have caused the disappearance of distinct ecotypes in all lakes except for Lake Superior (Agassiz, 1850; Brown et al., 1981; Goodier, 1981). There are currently four distinct ecotypes of lake charr still present in Lake Superior, leans, siscowet, humpers, and redfins (Muir et al., 2014). Ecotypes can be distinguished morphologically, and express different but often overlapping depth, habitat, and diet preferences. Studies that have examined the genetic isolation of lake charr ecotypes in Lake Superior using neutral markers demonstrate that genetic differentiation among ecotypes is often lower than spatial genetic differentiation (Baillie et al., 2016 et al., 2016; Page et al., 2004). Even across a small geographical range (approximately 10 km2), genetic variation among ecotypes appeared to be better described along a gradient of depth from shallow to deep, than by divergent ecotypes defined by morphometry (Baillie et al., 2018; Baillie, Muir, Hansen, et al., 2016). Furthermore, differentiation between ecotypes appears to have declined in recent years suggesting that there may be renewed gene flow among ecotypes (Baillie et al., 2016). The loss of genetic distinction could mean that the phenotypic differences between lake charr ecotypes in Lake Superior are also beginning to erode. However, it is difficult to make this conclusion without a better understanding of the relative role that genetics and environment play in determining key traits associated with each ecotype.

Many studies have attempted to dissect lake charr diversity to determine whether phenotypic variation is a result of hard‐wired genetic differences or phenotypic plasticity associated with environment. One of the most persistent phenotypes that differentiate lake charr ecotypes is muscle lipid content (Sitar et al., 2020). Ecological studies suggest that lipid content and other traits that differentiate lake charr ecotypes are related to depth habitat, foraging strategy, and buoyancy whereby ecotypes are specialized to different diet niches within the ecosystem (Muir et al., 2014). Many of the same traits and selective forces identified in Lake Superior, including buoyancy, are also important factors that distinguish sympatric lake charr ecotypes that have evolved in other lakes (Zimmerman et al., 2006, 2007). However, the exact ecological selection pressure is unknown. Genomic studies of wild populations using putatively non‐neutral loci, have found limited evidence of strong diversifying selection among ecotypes in lakes with two or more sympatric lake charr ecotypes, but do suggest that differences may be polygenic (Perreault‐Payette et al., 2017). Muscle lipid content, coloration, and body shape, all of which appear to differentiate among lake charr ecotypes to some degree (Muir et al., 2014), are complex traits probably derived from multiple genes (i.e., polygenic). Furthermore, Salmonines experienced a full genome duplication potentially increasing the number of loci which can potentially influence any given trait (Salisbury & Ruzzante, 2021). Therefore, allele frequency differences at loci linked with key traits may be difficult to distinguish from spatial genetic variation, thereby masking evidence of selection in studies of wild populations (Perreault‐Payette et al., 2017; Santure & Garant, 2018; Smith et al., 2020). Continued investigation of the genetic association of growth and lipid content in a common garden setting could help to identify possible targets of selection and further determine the heritability of lake charr ecotypes by controlling for environmental stochasticity.

Common garden rearing studies on lake charr ecotypes are beginning to provide important insight into the roll that genetics plays in determining phenotypic variation among and within lake charr ecotypes (McDermid et al., 2007, 2013; Smith et al., 2020). For example, truss measurements (measurements of head and body shape) that are used to differentiate wild lean and siscowet lake charr (Moore & Bronte, 2001) also differed between laboratory‐reared leans and siscowets (Goetz et al., 2010). Ecotypes also appear to fundamentally differ in the way they store and process energy as shown by differentiated glycogen, lipid and glucose levels of lake charr reared in captivity (Goetz et al., 2014). However, the most dramatic observation in the laboratory‐reared lake charr was that siscowets had significantly higher lipid in the muscle than leans in year one, and this difference increased further into years two and three (Goetz et al., 2010, 2014). Transcriptomic analysis indicated that gene groups involving lipid synthesis, metabolism, and transport were differentially expressed among lean and siscowet lake charr reared in a common garden environment, further supporting a role of lipids in determining lake charr ecotypes (Goetz et al., 2010). These studies proved that phenotypic and metabolic differences between siscowet and lean lake charr ecotypes were largely not mediated by environment, and therefore heritable. However, without data from an additional F1 generation or genetic data from nuclear DNA, these studies could not determine if phenotypic differences were the result only of transcriptional differences or maternal effects associated with the parental lineages.

The common garden rearing study first described in Goetz et al. (2010) was started in 2006 using offspring derived from wild parents and is now in its 15th year. The lean and siscowet lake charr used in this study exemplify the two depth niche extremes in Lake Superior; leans being considered a shallow water ecotype whereas siscowet are deep water ecotypes (Bronte et al., 2003; Sitar et al., 2008). The second generation of leans and siscowets are now nine years old and the objective of the present study was to test the hypotheses that (1) growth, development, and lipid content variation among full ecotype and half ecotype F1 crosses is stable through time and (2) phenotypic differences among ecotypes are genetically regulated by an overlapping set of genes and therefore heritable. Preliminary results of the F1 crosses indicated that lipid differences between leans and siscowets persisted into the first 3 years of the second generation (Goetz et al., 2014). Based on these results, we predict that growth and lipid differences among crosses would stabilize as individuals approached maturity if traits were fully heritable and not a consequence of maternal effects. Due to the complex nature of lipid and growth phenotypes, we predicted that we would find evidence of additive and polygenic variation present within both ecotypes and all four experimental crosses of F1 individuals.

MATERIALS AND METHODS

The lake charr used in the current investigation have been part of a common garden rearing study investigating the basis of the phenotypic differentiation of lean and siscowet ecotypes. The phenotypic differentiation of the parental lines (P1) was described in Goetz et al. (2010) and these lines were used to generate the F1 crosses described in the current study (Figures S1 and S2). The crossing and rearing scheme for producing the F1s is described in detail in Figure 1. Briefly, the original lean and siscowet P1 lines were derived from gametes of wild adult fish obtained in 2006 from Lake Superior. The fertilized eggs and subsequent juveniles and adults were reared under identical environmental conditions from 2006, initially at the Great Lakes WATER Institute (2006–2010: GLWI, School of Freshwater Sciences, University of Wisconsin – Milwaukee) and later at the Northern Aquaculture Demonstration Facility (2011–2019: NADF, University of Wisconsin – Stevens Point) where the F1 generation was produced and reared to adults. In the fall of 2011, at least 20 matings (two males × one female) of lean × lean [L × L], siscowet × siscowet [S × S], siscowet (female) × lean (male) [S × L], and lean (female) × siscowet (male) [L × S] were made. Eggs of individual crosses were kept separate in vertical Heath incubators and following swim‐up, approximately 2000 fry in total were randomly collected from at least 10 matings per cross and placed in four separate but identical tanks, each containing fish from a different cross and supplied with flow‐through water from the same source. Feed habituation to artificial trout diet (EXTR 450‐Rangen, Inc., Buhl, Idaho) was carried out in the tanks under identical conditions. Fish were progressively culled to 500 fish/tank and, after one year, 220 fish from each cross were randomly sampled from each tank and tagged with Passive Integrated Transponder tags. After tagging they were comingled and held in two identical tanks with flow‐through water and then temporarily in a raceway until they were 3 years old at which time they were distributed to eight identical tanks. At this same time, fin clips were collected for later genetic analysis. Each tank contained approximately 20 fish/cross (80 fish total/tank). They were reared in these tanks under the same environmental and feeding conditions until 2020. Natural mortality occurred approximately evenly across all crosses throughout the duration of the experiment at a low rate (~1.0% per year). In 2017, the number of tanks was reduced from 8 to 7 when fish from one tank were used for a separate analysis. All fish involved in these studies were reared and assessed under protocols approved by the University of Wisconsin‐Steven’s Point and the University of Wisconsin‐Milwaukee Animal Care and Use Committees.

FIGURE 1

Crossing and rearing scheme to produce the F1 lake charr used in the current study. F1 lake charr lines were produced from lean and siscowet parents that were produced originally in 2006 using wild gametes obtained from Lake Superior (Goetz et al., 2010). In the fall of 2011, at least 20 single pair matings (two male ×one female) of lean ×lean [L × L], siscowet ×siscowet [S × S], siscowet (female) × lean (male) [S × L], and lean (female) × siscowet (male) [L × S] were made. Eggs from these matings were incubated and the fry were subsequently reared to adults at the Northern Aquaculture Demonstration Facility (University of Wisconsin–Stevens Point). Eggs of individual crosses were kept in separate trays in vertical Heath incubators supplied with flow‐through water at 8–10°C. Following swim‐up, approximately 2000 fry in total were randomly collected from at least 10 matings per cross and placed in four separate but identical tanks (circular 946 L), each containing fish from a different cross and supplied with flow‐through water (7.6–8°C) from the same source. Feed habituation to artificial trout diet (EXTR 450‐Rangen, Inc., Buhl, Idaho) was carried out in the tanks under identical conditions. Fish were progressively culled down to 500 fish/tank and, after one year, 220 fish from each cross were randomly sampled from each tank and tagged with passive integrated transponder tags. After tagging they were comingled and held in two 1960 L tanks with flow‐through water (7.6–8°C) and then temporarily in a raceway (7.6–8°C) until they were 3 years old at which time they were distributed to 8–1960 L circular tanks. In 2017 the number of tanks was reduced to 7. Each tank contained equal numbers of all four crosses. They were reared in these tanks under the same environmental and feeding conditions (7.9–8.5°C) until 2020

The influence of holding tank on phenotype was evaluated following the end of the experiment in 2019 using ANOVA with tank, year, and cross as explanatory variables (Figures S3–S5). There were no significant impacts of tank on lipid content (p = .12), weight (p = .09) and a very minor impact tank on length (p = .05) which was not significant following Tukey’s HSD corrections. Therefore, we did not consider tank as a factor in the rest of analyses. Furthermore, while fish were assessed for sex, sex was not considered as a significant covariate as numerous studies have reported that lake charr are not sexually dimorphic (Gunn, 1995; Martin, 1980; Muir et al., 2014).

Length, weight, and lipid assessment

In the fall (September–October) from 2012–2019, all fish within the F1 generation crosses were assayed for length (nearest 1.0 mm) and weight (nearest 1.0 gm [years 2–3] or 10 gm [years 4–9]). From 2015–2019, somatic lipid content (as a percent) was measured for each fish using a battery powered, handheld microwave oscillator (Distell Model 692 Fish Fatmeter, Distell Inc.). The fatmeter emits a low‐powered microwave (2 GHz, 2000 MHz, power 2 mW) that interacts with water within the somatic tissues and uses the inverse relationship between water and lipid to estimate the lipid concentration (as a percent) in the tissue (Crossin & Hinch, 2005; Kent, 1990). Fatmeter readings taken on the Research‐1 setting were collected on each fish at a site that was in the epaxial muscle mass just posterior to the head (site “S1” as shown in Sitar et al., 2020). We have previously shown that fatmeter readings and biochemical estimates of lipid content in lake charr muscle are highly correlated (Sitar et al., 2020), indicating that the fatmeter readings provide a reliable index of muscle lipid levels in lake charr. We also observed high correlations between the fatmeter readings behind the head (used in the current study) and the readings from other sites on the body (Sitar et al., 2020), indicating that the site just behind the head can be used as a reliable proxy for lipid at other muscle sites. From 2012 to 2014 the fish were too small to use the fatmeter on, but a subsample (n2012 = 12; n2013‐2014 = 20) of fish from each cross were lethally sampled and a muscle sample (cross‐section of the epaxial muscle behind the head and anterior to the dorsal fin) was excised and used for lipid analysis by Soxhlett extraction as described in Goetz et al. (2014) and Sitar et al. (2020). At the conclusion of the study, a single lipid content and Fultons condition factor (Equation 1) score was calculated for each individual as a simple average of all years of available data. This value was used as the phenotypic measure for genetic association tests.

Differences in lipid content, condition, and growth

Variance in lipid content based on biochemical analysis (years 2012 to 2014) and fatmeter readings (years 2015 to 2019), length, weight, and Von Bertlanffy growth parameters (Equation 2) among crosses was evaluated in R using custom scripts implementing the FSA and nlstools R packages (Baty et al., 2015; Cooper & Carlander, 1951; Ogle, 2014; R Core Team, 2019). Analysis of variance (two‐way ANOVA) was used to test for differences in lipid content, length, and weight among crosses with cross and year as the principal factors following qqplots confirming assumptions of normality (Figures S6–S9). To account for the change in lipid measurement from biochemical analysis to fatmeter readings the first three years of data were analysed separately from the last five using identical custom R scripts. Following two‐way ANOVA analysis, Tukey’s honestly significant tests were used to identify intra‐annual differences among crosses. Variation in growth among crosses was evaluated by comparing Von Bertalanffy growth parameters (L∞, K, and t0; Murphy & Willis, 1996). The Von Bertalanffy growth function is used to estimate average growth based on size‐at‐age data which produces three parameters, L∞, which is the asymptotic average length, K which is the Brody growth coefficient (units per year), and t0, which is a modeling artefact representing age at length 0. Parameters were estimated independently for each cross using the FSA package in R (Ogle, 2014). Starting parameters were identified using the vbStarts function and the typical growth model parameter. Typical growth models were then fit using nls function and bootstrapped with 999 bootstrapping iterations. Differences in growth parameters were inferred based on the overlap in bootstrapped confidence intervals.

Parentage

Parentage of all individuals was assigned using genotypes from a six‐microsatellite marker panel amplified and genotyped for all parents and offspring using previously designed primers for Loci SnaMSU 01, 03, 06, 10, 11, and 12 (Rollins et al., 2009) and PCR and genotyping protocols outlined in Rosauer et al. (2011). Pedigrees were reconstructed using Colony2 Version 2.0.6.5 and a full‐likelihood approach specifying maternal polygamy without inbreeding and no sibship prior (Jones & Wang, 2010). Marker error rate and probability of null alleles was assumed to be 0.0001 and genotypes for all putative dams and sires were provided and the probability that a parent was in the data set was set to 1. Parentage for each cross was run separately using identical parameters and parents were defined as the highest probability parent pair. Suggested pairings were then filtered to include only pairings of individuals known to have been spawned together.

DNA sequencing

Restriction site associated DNA sequencing (RAD‐seq) was conducted on 74 parents and 542 F1 offspring using SbfI and bestRAD library protocols outlined in Ali et al. (2016). An initial library was sequenced at BGI America on one lane of a HiSeq4000 and the remainder of the libraries were sent to Novogene where they were sequenced on seven HiSeq4000 lanes for paired‐end 150 sequencing. At the time of data collection, no lake charr genome assembly was available, therefore identification of SNPs and genotyping were conducted in STACKS v.2.3 using the de novo assembly pipeline (Rochette et al., 2019). Samples were demultiplexed with process_radtags (flags = c, −q, −r, −t 140). Stacks of similar sequences (loci) for each individual were identified with ustacks (flags = −m 3, −M 5, −H –max_locus_stacks 4, –model_type bounded, –bound_high 0.05) and a catalogue of putative loci was generated based on sequences from the parents. Individual stacks were then aligned to the catalogue in sstacks, and genotypes for all putative SNPs were assigned using gstacks. Finally, a datafile containing genotypes for all SNPs with a minor allele count greater than two and all individuals was generated using populations and subsequently filtered in VCFTools (Danecek et al., 2011). To ensure that paralogues did not influence our findings, we ran HDPlot on unfiltered data for all crosses and removed loci identified as potential paralogues (McKinney et al., 2017). Parameters for this analysis were set by visually choosing threshold values for read depth ratio and proportion of heterozygotes that identified the loci conforming to theoretical expectations for singletons (McKinney et al., 2017). Once putative paralogues were removed, quality filtering was conducted in VCFTools hierarchically based primarily on SNP genotyping rate. Because the individuals sampled are from an experimental population, we expected that loci with accurate genotype calls may be out of Hardy‐Weinberg Equilibrium and have elevated linkage disequilibrium and therefore did not use these metrics to filter loci. First, all individuals with more than 80% missing data were removed from the data set. Second, SNPs missing more than 10% of genotypes were removed. Finally, for loci that contained more than one SNP, the SNP with the highest minor‐allele‐frequency was retained leaving a single SNP per‐locus. All sites that passed the above thresholds were retained in downstream analysis. All bioinformatic analysis was conducted using the Turing High Performance Computing cluster at Old Dominion University, Virginia.

Genome‐wide cross differences

Genome wide association studies (GWAS) are sensitive to population structure, and relatedness (Santure & Garant, 2018). When our experiment began in 2006, siscowet broodstock were collected from a different spawning shoal approximately 320 km east from where experimental lean lake charr were collected (Figure S1). Therefore, it is likely that neutral population structure not linked to ecotype could bias our analysis (Baillie et al., 2018; Baillie, Muir, Hansen, et al., 2016; Perreault‐Payette et al., 2017). To assess and visualize total genetic variation among crosses, we conducted principal component analysis (PCA). Missing genotypes were filled by imputing the most common genotype at the locus within each cross and principal component analysis was run on the complete data set of SNPs in the Vegan R package (Oksanen et al., 2020).

Identification of trait associated candidate loci and QTL analysis

To identify loci associated with phenotype we first conducted a genome wide association test (GWAS) using PLINK v.1.9c (Purcell et al., 2007). To account for relatedness, PLINK was conducted using the family‐based association analysis for quantitative traits module and only individuals with known ancestry and sex. This approach uses linear regression of phenotype on genotype and permutation tests to account for family structure. Individuals from all crosses were analysed together but a strictly within‐family design was used to calculated permuted‐p‐values. Significance of associations was assessed using both Benjamini and Hochberg (1995) adjusted p‐values and Bonferroni corrected p‐values. Following initial assessment, conclusions were identical following both corrections and so the more conservative Bonferroni correction is reported.

To account for neutral population structure, following initial GWAS in PLINK, we focused our genetic association analysis on each cross separately to identify SNPs that explained a disproportionate amount of variation in lipid content and condition factor among families within each cross. Putative genetic association with lipid content and condition factor within each cross was evaluated using redundancy analysis (RDA). Prior to modeling, potential explanatory variables (family, mean lipid content, weight, and length) were evaluated by assessing pairwise‐correlations using pairs.panels() in the pych R package (Revelle, 2021) and by calculating variance inflation factors (VIF) to limit correlations among predictors that can lead to overfitting of the model. Based on these preliminary analyses, length and weight were combined into a single variable “condition” as they were highly correlated (R2 = 0.93). After consideration, family was also dropped as an explanatory variable. For many families, VIF scores were high (>10) which seemed to produce overfit models and spurious correlations to family, and not phenotypes. Furthermore, family assignment is completely dependent on individual genotypes and therefore breaks the independence assumption central to RDA (Xia, 2020). Therefore, the final models included mean lipid content and condition as the only two explanatory variables which constrained the ordination to identify only variation explained by phenotype, regardless of family. Each overall RDA model as well as each RDA axis was assessed for significance using ANOVA and 999 permutations in Vegan. Suggestive candidates were classified as SNPs with loading scores greater than three standard deviations from the mean SNP loading score for all RD axes (Forester et al., 2018). Suggestive candidate SNPs were then screened further to retain only SNPs that were either positively or negatively correlated with either lipid content or condition factor in all four crosses. Genes closest to each candidate SNP were obtained by aligning candidate sequences to the lake charr genome using BWA and identifying the nearest gene with the bedtools v.2.22.0. closest tool (Quinlan & Hall, 2010; Smith et al., 2021). Based on Smith et al. (2020), genes within 50 kb of candidate loci were defined as potential targets of observed trait variation and researched using UniProt (version 10 April 2018) and QuickGO (version 6 June 2021) and a subsequent literature search. To evaluate the additive effect of candidates on phenotype, we identified which allele at each of the candidate SNPs was positively correlated with phenotype and then calculated the sum of positive effect alleles (PEAs) that each individual had for lipid content and condition separately to create a polygenic score (Chatterjee et al., 2016). The relationship between the sum of PEAs and phenotype was evaluated using a Pearson correlation conducted separately for each cross and the overall influence of PEAs and estimates of explained variance were determined using a 2‐way ANOVA (linear equation: mean phenotype ~sumPEA + cross). Approximate normality of variables was assessed using qqplots prior to analysis. This approach assumes a model of incomplete dominance which is probably not the case for all markers. However, the presence of a positive correlation was considered to be additional support that loci are associated with phenotype since if SNPs were false positives, there should be no additive effect.

Because individual lipid content and condition are quantitative traits influenced by many loci, we also mapped loci and traits to a lake charr linkage map and scanned the genome for potential quantitative‐trait‐loci (QTL; Smith et al., 2020). The position of each SNP in cM was estimated by aligning sequences for all RAD loci in our data set and sequences of all QTL loci reported by Smith et al. (2020) in Table S1 to the lake charr genome with BWA (NCBI: GCF_016433055.1; Smith et al., 2021). All alignments with a score less than 60 were removed, as well as any locus that aligned to a different chromosome than reported in Smith et al. (2020). The position in cM along the female linkage map was estimated for our SNPs by using the predict function in R and an equation for a loess curve fitted from the Smith et al. (2020) alignment and linkage map SNP position for each linkage group. Once the position of SNPs was estimated, QTL mapping was conducted following almost identical procedures as outlined in Smith et al. (2020) with the qtl2 R‐package (Broman et al., 2019) and individuals from all crosses with known pedigrees and a family size greater than five individuals. Pseudomarkers were added to the map using insert_pseudomarkers with a step size of 1 cM and calc_genoprob was used to calculate genotype probabilities. Markers were thinned with calc_grid and probs_to_grid and used to create a kinship grid with calc_kinship which was used to account for confounding effects of family. Scans for QTL were conducted using scan1 and then find_peaks were used to identify suggestive peaks (threshold = 3, peakdrop = 3, drop = 2) and to estimate 95% credible intervals around each peak (prob = 0.95, peakdrop = 2, threshold = 3). Significance LOD (log likelihood of odds) thresholds (alpha = .05) were generated for each trait by permuting a null distribution 999 times using scan1perm. All genes within the 95% CI for significant QTL were identified as potential targets of trait variation within lake charr. All figures and data summaries were conducted in R using ggplot2, ggpubr, tidyverse, and ggsci packages (Kassambara, 2020; Wickham, 2009; Wickham et al., 2019; Xiao, 2018). A summary of essential code used during analysis can be found on GitHub (https://github.com/peuclide).

RESULTS

Phenotype differences

Variance of all phenotypic measurements (lipid content years 2012 to 2014, lipid content years 2015 to 2019, weight and total length) significantly diffed among crosses and years and showed a significant cross by year interaction (Table 1). Full siscowet crosses were the most lipid rich (2012 to 2014 average = 44.5%; 2015 to 2019 average = 55.4%) and full lean crosses were the least lipid rich of the four crosses (2012 to 2014 average = 22.5%; 2015 to 2019 average = 37.0%; Figure 2). Both Lean (Female) × Siscowet (Male) and Siscowet (Female) × Lean (Male) half crosses had lipid content that was intermediate to full crosses (2012 to 2014 average = 33.8% and 31.8%; 2015 to 2019 average = 47.8% and 48.5%) and did not significantly differ in most years. While there was a significant cross by year interaction, Tukey’s HSD tests indicated that differences in lipid content among crosses were mostly stable from one year to the next (Figure 2). The pattern in lipid content only differed in the first (2012) and last (2019) years of the experiment. In both years, Siscowet (Female) × Lean (Male) lake charr appeared to have higher lipid content than Lean (Female) × Siscowet (Male) lake charr, however this effect was only significant in 2019 according to Tukey’s tests (2012 ptuk = 0.98; 2019 ptuk < 0.01). This is probably due to the small sample size used for biochemical analysis in 2012 (Table 1).

TABLE 1

Results from two‐way analysis of variance comparing lipid content, weight, and length among experimental F1 lake charr crosses reared in captivity for eight years. Biochemical lipid content readings (% lipid content) were used from 2012 to 2014 followed by fatmeter readings (% lipid content) starting in 2015, and therefore data reported separately

Effect	df	Residuals	F	p‐value
Lipid content (2012–2014)
Cross	3	196	210.4	5.1 × 10−61
Year	2	196	52.6	5.1 × 10−19
Cross:year	6	196	11.3	7.5 × 10−11
Lipid content (2015–2019)
Cross	3	2,362	539.0	<1.0 × 10−100
Year	4	2,362	50.7	5.7 × 10−41
Cross:year	12	2,362	3.4	6.0 × 10−5
Length (mm)
Cross	3	4,261	110	2.0 × 10−68
Year	7	4,261	10387.0	<1.0 × 10−100
Cross:year	21	4,261	14.4	1.0 × 10−49
Weight (g)
Cross	3	4,262	56.8	5.3 × 10−36
Year	7	4,262	3018.0	<1.0 × 10−100
Cross:year	21	4,262	15.1	2.0 × 10−52

FIGURE 2

Variation in percent lipid content of muscle tissue (a) length in millimeters (b) and weight in grams (c) for F1 lake charr crosses reared under identical conditions for nine years. Percent lipid content measured in the years 2012, 2013, and 2014 (left of the dotted line) were estimated biochemically while percent lipid content for the remaining five years (2015–2019) was measured using a fatmeter reading. Percent lipid content measured biochemically was previously found to be tightly correlated with fatmeter readings Sitar et al. (2020) and therefore results are plotted together to simplify interpretation. Results of Tukey’s HSD tests are shown as letters at the top of each boxplot with different letters indicating significant intra‐annual difference between crosses. Statistical differences in lipid content among crosses for years 2012 to 2014 and 2015 to 2019 were analysed separately to reduce bias associated with the change in methodology. Lipid content, length, and weight also differed significantly by year (Table 1)

Length and weight cross differences followed almost identical patterns of each other throughout the duration of the experiment (Figure 2). Slight differences in length and weight appeared at age 2 (2013) and primarily occurred between full lean and full siscowet crosses. Length and weight of half crosses were most often intermediate of full crosses; however, at age 2 (2013) and 3 (2014) Lean (Female) × Siscowet (Male) crosses appeared to be larger than the other crosses. Patterns in both weight and length began to stabilize around age 4 (2015). Starting at this age, half crosses largely reflected their maternal lineage whereby Lean (Female) × Siscowet (Male) individuals were similar in length and weight to full lean crosses and Siscowet (Female) × Lean (Male) individuals were similar in length and weight to full siscowet crosses. However, this pattern was not statistically supported in all years (Figure 2).

Full siscowet crosses grew more quickly (KTL) than all other crosses in the first three years but also reached maximum size (L∞) more quickly (Figure 3). Von Bertlanffy growth curves reflected the same patterns as lipid content and condition whereby half crosses had intermediate growth of individuals from the two full crosses. Full siscowet crosses had the smallest L∞ (mean = 548, 95% CI: 542–556 mm) followed closely by Siscowet (Female) × Lean (Male) individuals (mean = 567, 95% CI: 560–574 mm), and Lean (Female) × Siscowet (Male) individuals (mean = 581, 95% CI: 575–588 mm) and then full lean crosses had the largest L∞ (mean = 624, 95% CI: 614–634 mm). Growth rate (KTL) was slowest for full lean crosses (mean = 0.38, 95% CI: 0.36–0.40), followed by Siscowet (Female) × Lean (Male) individuals (mean = 0.42, 95% CI: 0.40–0.43), and Lean (Female) × Siscowet (Male) individuals (mean = 0.48, 95% CI: 0.46–0.50) and then full siscowet crosses had the highest growth rate (mean = 0.53, 95% CI: 0.50–0.55). Estimated t0 ranged from 0.33 to 0.50 and 95% confidence intervals overlapped for all crosses (Figure 3d).

FIGURE 3

Estimated growth metrics for F1 lake charr crosses reared under identical conditions for nine years. Predicted length at age is summarized in (a) whereby each line represents the predicted growth trajectory of each cross and the width of each line corresponds to bootstrapped lower and upper confidence intervals. (a–d) Summarize the variance of estimated von Bertalanffy growth parameters: projected maximal growth (L∞), growth rate (K), and size at age‐0 (t0)

Pedigree

Parents were successfully assigned in Colony to 521 of the 542 individuals. Across the four crosses a total of 45 families from 24 dams and 45 sires. Family success was variable with nine families producing 52% of all progenies. Each cross contained a total of 10–14 families, but two to four families produced greater than 50% of offspring in all crosses (L × L = 3 families; L × S = 2 families; S × L = 2 families; S × S = 4 families).

RAD‐sequencing

RAD‐sequencing produced 3,584,739,804 reads of which 3,484,674,503 were retained resulting in an average effective coverage per‐individual of 31.9 (SD = 21.3). A total of 347,213 putative SNPs with a minor allele count greater than two were identified. Prior to SNP filtering, eight individuals were removed that had a genotype rate less than 80%. Putative paralogues were identified and removed based on z score greater than 5 or less than –5 and or with a heterozygosity greater than 0.5 (Figure S10). Next, all SNPs that genotyped in less than 90% of individuals were removed. Finally, to limit the influence of physical linkage among SNPs on the same RAD‐tag, only the SNP with the highest minor allele frequency at each locus was retained. Following all filters, a total of 19,714 SNPs were retained for analysis.

Phenotype heritability

Genotypic variance was primarily explained by cross (Figure 4a). The variance explained by PC axis 1 was 6.4% and strongly isolated S × S from L × L crosses while S × L and L × S individuals showed intermediate values for PC1 as would be expected of a F1 hybrid. PC axis 2 explained 4.8% of the variance and appeared to differentiate S × L and L × S families. Additional PC axes appeared to further delineate among families. No SNPs were significantly associated with either lipid content or condition when all crosses were analysed together in PLINK, and no chromosomes appeared to have elevated SNP associations following permutations to account for family (Figure S11). This analysis indicated that there were no clear patterns of high‐effect loci linked to either phenotype, and we therefore chose to conduct the rest of our analysis within each cross using a different approach (RDA) to identify more subtle patterns of allele frequency differences within crosses while accounting for neutral structure. Following initial tests for model assumptions outlined in the methods (Figure S12, S13), we conducted RDAs using the complete data set of SNPs and lipid content, with condition factor as variables. The RDA models for each cross was significant (p < .01). When tested separately, the first RD axis was significant for all crosses (p <.05) while the second RD axis was not (L × L p = .14; L × S p < .01; S × L p = 0.15; S × S p = .09; Table S14). Approximately 8% of all SNPs were identified as suggestive candidates in each cross (L × L = 1385 SNPs, L × S = 1379 SNPs, S × L = 1430 SNPs, S × S = 1435 SNPs). However, only 51 SNPs (0.25%) of the total 19,714 SNPs had suggestively high loading scores in all four crosses (Figure 4b). Of these 51 SNPs, 13 candidates contained alleles that were consistently positively or negatively correlated with lipid content and three candidates contained alleles that were consistently positively or negatively correlated with condition factor (Table 2). One SNP (CLocus_28583) was a candidate for both lipid content and condition factor. All 15 candidates were successfully aligned to the lake charr genome. Of these, 14 aligned to lake charr chromosomes and the remaining one SNP aligned to an un‐scaffold read (Table 2).

FIGURE 4

Genetic variation and putative phenotype associations identified using redundancy analysis (RDA) based on 19,714 single nucleotide polymorphisms (SNPs) identified with RAD‐sequencing of experimental F1 lake charr crosses (legend; F, female; M, male). (a) Genetic variation among crosses was summarized with principal component analysis of the first two eigen vectors and all SNPs. (b) Venn diagram of SNPs that explained a disproportionate amount of variance in lipid content or condition based on having an RDA loading score greater than 3 standard deviations from the mean. The 51 SNPs that were outliers in all crosses were further interrogated and only SNPs containing alleles that were positively or negatively correlated with phenotype in all crosses were retained as candidates for phenotype association. (c) and (d) show the relationship between then number of positive effect alleles (PEAs) at these candidates (maximum number of lipid content PEAs = 26, maximum number of condition PEAs = 6) and individual lipid content (c) and condition (d). Linear models are shown as lines with 95% confidence intervals of the relationship is shown as grey boxes around each line

TABLE 2

Candidate loci that were found to explain a disproportionate amount of variance in lipid content and condition based on redundancy analysis (RDA) that also contained alleles with a positive relationship with either lipid content or condition in all four crosses. Candidates were aligned to a draft lake charr genome and linkage map and annotated based on the closest gene or predicted gene to the alignment position in the lake charr genome

Trait	Locus	Chromosome	Position (Mb)	Estimated position (cM)	Closest gene	Distance from gene (Kb)	Gene type
Lipid	CLocus_69893	Chr1	54.17	37.65	LOC120045774	0	Diacylglycerol kinase zeta‐like isoform X5
Lipid	CLocus_43554	Chr5	61.22	6.08	LOC120047228	2.1	Meprin A subunit beta‐like
Lipid	CLocus_42626	Chr5	18.28	95.99	LOC120048021	0	LIM and calponin homology domains‐containing protein 1‐like
Lipid	CLocus_14988	Chr13	22.64	45.75	LOC120058335	0	Integrin beta−5‐like
Lipid	CLocus_91657	Chr21	19.14	44.47	LOC120066156	0	SRSF protein kinase 2‐like
Lipid	CLocus_23471	Chr22	7.36	35.52	LOC120017377	3.2	Pyruvate dehydrogenase E1 component subunit alpha, somatic form, mitochondrial‐like
Lipid	CLocus_47288	Chr22	7.52	35.95	ptpn20	18.9	Tyrosine‐protein phosphatase non‐receptor type 13
Lipid	CLocus_67283	Chr24	12.38	47.19	LOC120019286	0	Solute carrier family 35 member F1‐like
Lipid	CLocus_33926	Chr26	8.53	26.79	fibcd1b	20.1	Fibrinogen C domain‐containing protein 1
Lipid	CLocus_97377	Chr35	28.04	35.17	LOC120029674	15.8	Ankyrin repeat domain‐containing protein 1‐like
Lipid	CLocus_91695	Chr38	15.56	6.37	LOC120031864	0	Protein mono‐ADP‐ribosyltransferase PARP12‐like
Lipid	CLocus_28583	Chr34	35.33	47.82	LOC120028506	211.4	Putative uncharacterized protein DDB_G0286901
Lipid	CLocus_73675	Chr34	37.43	48.29	LOC120028514	54.6	Protocadherin Fat 3‐like
Condition	CLocus_48397	Chr23	28.54	9.27	trim59	0	Tripartite motif‐containing protein 59
Condition	CLocus_50919	NW_024058375.1	0.18	–	–	–	–
Condition	CLocus_28583	Chr34	35.33	47.82	LOC120028506	211.4	Putative uncharacterized protein DDB_G0286901

All crosses contained individuals with at least three and a maximum of 22 PEAs and qqplots indicated that the distribution of PEAs in each cross were approximately normal, though the distribution was slightly left skewed for full the lean cross individuals. Pearson’s correlations indicated that there was a positive relationship between the number of PEAs and both muscle lipid content and condition in each of the crosses (Figure 4c,d; Table 3). Based on R2 values, PEAs explained 5 to 19% of the variance in lipid content within each cross and 3 to 5% of the variance in condition (Table 3). When comparing across all crosses (i.e., a two‐way ANOVA with PEAs and cross as factors), PEAs explained very little of the variance in lipid content and condition. The effect of PEAs on overall lipid content was not significant (F1,529 = 0.3; p = .5), and explained only 2% of the total variance in muscle lipid content while experimental cross explained 62% of the variance (F3,529 = 287.6; p < .001). Condition was also significantly influenced by the number of PEAs present (F1,529 = 22.7; p < .001) but also primarily differed by cross (F3,529 = 27.0; p < .001). However, when the amount of variance in condition explained by PEAs and cross was compared to the amount of variance in lipid content explained by PEAs and cross, the effect was modest and PEAs explained 4% of the variance while cross explained 13%.

TABLE 3

Pearson’s correlation summary statistics between the sum of positive effect alleles calculated for each individual and individual phenotype shown in Figure 3c,d

Cross	Trait	Estimate	Lower confidence interval	Upper confidence interval	Test statistic	Degrees of freedom	p‐value	R2
L × L	Lipid	0.23	0.06	0.39	2.65	123	.009	0.05
L × S	Lipid	0.30	0.14	0.44	3.76	141	.000	0.09
S × L	Lipid	0.43	0.28	0.56	5.56	135	.000	0.18
S × S	Lipid	0.43	0.28	0.56	5.41	130	.000	0.18
L × L	Condition	0.22	0.05	0.38	2.52	123	.013	0.05
L × S	Condition	0.18	0.01	0.33	2.12	141	.035	0.03
S × L	Condition	0.22	0.05	0.37	2.58	135	.011	0.05
S × S	Condition	0.21	0.04	0.36	2.41	130	.017	0.04

Candidates SNPs spanned 12 chromosomes with Ch5, Ch22, and Ch34 each containing two candidates. Seven of the candidates intersected either known or predicted genes and five were within 50 kb of a known or predicted genes. The genes that contained candidate loci were broadly associated with binding of metal ions and cell cycling through apoptosis and growth. CLocus_69893 aligned to a gene isoform which has a broad range of functions associated with cellular function and cell signalling (Van Blitterswijk & Houssa, 2000). CLocus_42626 aligned to a gene isoform which is involved in calcium ion binding (GO:0046872). CLocus_14988 aligned to a gene orthologue that plays a role in cell signaling pathways and has been found to be important for angiogenesis, vasculogenesis, hematopoiesis, and bone formation in mice (Gat et al., 2001). CLocus_91657 aligned to protein isoform which is important of mediating neuronal cell cycles and development (Bustos et al., 2020). CLocus_67283 aligned to a protein which is a transmembrane transporter protein involved in endoplasmic reticulum homeostasis (Nishimura et al., 2009). CLocus_91695 aligned to a protein which is another metal ion binding protein associated with basic cell functions including DNA repair and apoptosis as well as expression of growth factor beta1 (MacPherson et al., 2013; Shao et al., 2018). This locus was also closest to the significant QTL identified on Chromsome 38. CLocus_48397 aligned to the trim59 protein which is another metal ion binding protein that is believed to be involved numerous cancers and immune response (McNab et al., 2011; Valiyeva et al., 2011).

One significant QTL peak was identified for lipid content on Chr38 (5.3 cM, 95% CI: 4.3–6.3, LOD = 5.7, p = .03). Furthermore, the upper confidence interval for the significant QTL peak on Chr38 was close to one of the 51 RDA candidate (CLocus_91695; 6.37 cM). The region of Chr38 located between the lower 95% CI of the significant QTL and CLocus_91695 contained 220 annotated genes and 58 SNPs identified by the present study. Several other suggestive peaks with LOD scores >3 were identified (Figure 5a). A large, suggestive QTL was identified on Chr1 (17.2 cM, 95% CI: 16.2–28.2, LOD = 5.1, p = .051) which was 9.4 cM away from candidate locus CLocus_69893 (37.6 cM). A smaller suggestive QTL on LG 5 (57.2 cM, 95% CI: 3.1–103.1 cM, LOD = 3.18, p > .9) included candidate locus CLocus_43554 (6.08 cM) and CLocus_42626 (96.5 cM). No statistically significant QTL were identified for condition, there were two large suggestive peaks that did not occur near RDA candidate loci and therefore were not further investigated.

FIGURE 5

Logarithm of odds ratios (LOD) values by estimated single nucleotide polymorphisms (SNP) position (cM) on the female lake charr linkage map reported in Smith et al. (2020) for quantitative‐trait‐loci scans for lipid content (a) and condition (b). The dashed red line corresponds to significance threshold (p < .05) and the solid purple line corresponds to the suggestive threshold of 3 used by Smith et al. (2020) to identify peaks putatively associated with traits. Open circles indicate predicted position in cM of candidate SNPs identified with redundancy analysis and described in Table 2

DISCUSSION

Taken together our results indicate that phenotypic differences in muscle lipid content between lake charr ecotypes were additive, polygenic, and heritable into the second generation of our common garden experiment. Key metabolic and growth differences between F1 siscowet and lean ecotypes could be observed within the first year and remained discernable and almost unchanged for the eight‐year experiment strongly indicating that lake charr ecotype is constrained by ancestry. Intermediate phenotypes of hybrid F1 individuals is consistent with additive trait differences and suggests that distinct ecotypes may be maintained through reproductive isolation. Furthermore, this reproductive isolation may be facilitated through divergent selection on lake charr morphology and metabolism. Our results also suggest that much of the genetic basis for lipid content may be polygenic. Other growth factors that differed within and among crosses, such as length and weight may also be partially genetically regulated, however the effect we observed was much smaller. While most of the variation in muscle lipid content and condition between ecotypes could not be explained by polygenic scores of candidate loci identified in the present study, we are able to conclude that there is a genetic basis for key traits that differentiate lake charr ecotypes.

Depth related trait differences

Sympatric phenotypic variation and resource polymorphisms have been observed numerous times among freshwater and marine organisms and are often associated with depth and foraging strategy (Farré et al., 2016; Moura et al., 2015). Charrs in the family Salmonidae exemplify this type of variability and frequently show a high degree of adaptability and phenotypic variation (Jonsson & Jonsson, 2001; Muir et al., 2016). However, for studies that have addressed the genetic and environmental basis for phenotypic variation in Arctic char (Salvelinus alpinus), some find predominantly genetic control (Skúlason et al., 1996) and others find predominantly environmental control (Adams & Huntingford, 2004). Lake charr show similar levels of adaptability and phenotypic variation as Arctic charr and sympatric ecotypes have evolved in several lakes that can be partitioned by depth strata (Chavarie et al., 2021). The morphology and prevalence of lean and siscowet lake charr ecotypes used in our experiment also diverge based on depth (Bronte et al., 2003; Sitar et al., 2008). Depth appears to describe a substantial amount of lake charr genetic diversity as well, suggesting that spatial ecology plays an important role in the maintenance of ecotype diversity (Baillie et al., 2018; Baillie, Muir, Hansen, et al., 2016). In line with what Skulason et al. (1996) found in Arctic Charr, our data suggest that key traits associated with different ecotypes, such as lipid content, is heritable, and supports the hypothesis that there is genetic regulation of lake charr phenotypes observed in the wild.

Depth and diet have been found to be the primary factor partitioning sympatric ecotypes of a freshwater fish multiple times (e.g. Dalziel et al., 2015; Rogers et al., 2002), but lipid content has rarely been discussed as an important trait. However, lipid content is one of the key traits that differentiates lake charr ecotypes in the wild (Eschmeyer & Phillips, 1965; Sitar et al., 2020) and one that persisted in common garden conditions in our study. Muscle lipid content variance was first described in the parental lineage (Goetz et al., 2010), but the variance remained substantial in the F1 generation presented here, whereby siscowets had 20% higher muscle lipid content than leans. Lipid levels vary greatly in fish (e.g. Devadason et al., 2016; Wang et al., 1990). Tissue lipid content is important for many reasons including buoyancy regulation (Clarke et al., 1984; Henderson & Anderson, 2002; Lee et al., 1975; Phleger, 1998), energy storage (Goetz et al., 2014; Sheridan, 1988, 1994), and for metabolism (Tocher, 2003). Given the strong pattern of heritability we observed, we hypothesize that lipid content serves an important function for lake charr and may be under selection.

Buoyancy regulation has been one of the main ecological mechanisms hypothesized to be driving lipid differences among lake charr ecotypes. In most teleost fishes, swim bladders allow fish to maintain position at different depths by moving gas in or out of the swim bladder when moving vertically under changing hydrostatic pressures (Pelster, 2021). Increased lipid in deep‐water lake charr ecotypes may make it easier for siscowet to regulate buoyancy without extensive gas exchange while undergoing large, and rapid changes of depth that are not frequently exhibited by lean lake charr (Eshenroder & Burnham‐Curtis, 1999; Henderson & Anderson, 2002). These vertical movements are probably involved in foraging and are extreme requiring fish to move hundreds of meters to the surface and then back again frequently within periods of 15–30 min (Binder et al., 2021; Jasonowicz et al., 2022; Stockwell et al., 2010). Repeatedly inflating and deflating a swim bladder during vertical movements would be metabolically taxing, and swim bladders may go unused when siscowet are in their primary habitat on the lake bottom (Goetz personal observation; McCune & Carlson, 2004). In contrast, lean lake charr spend much more of their time in pelagic waters where swim bladders would be beneficial for short vertical movements and for maintaining a relatively stationary depth (Pelster, 2021). The heritability of lipid content may therefore be a consequence of divergent selection acting on the different foraging and behavior needs of lean and siscowet lake charr.

Depth gradients may also influence growth and body size differences observed between lake charr ecotypes in our study. In the wild, siscowet lake charr reached a smaller asymptotic length and contained smaller Brody growth coefficient (K) than lean lake charr (Hansen et al., 2016). However, this is somewhat different from what we found. F1 siscowet in the current investigation did express smaller asymptotic growth, but our estimate of the Brody growth coefficient was reversed, and higher for siscowets than leans (i.e., siscowet may have grown more rapidly in captivity than they did in the wild or lean grew more slowly). This pattern was the same for the parental lineage, and therefore probably inherited (Goetz et al., 2010). One hypothesis is that constant experimental temperature conditions (7.5–8°C year‐round) led to faster or slower growth of one or both ecotypes. Thermal habitats for wild siscowets have not been directly recorded but would have to be 3–5°C year‐round when they were on or near the bottom based on the temperature profiles of Lake Superior (Titze & Austin, 2014). In contrast, leans occupy temperatures between 8 and 12°C during the summer closer to the surface (Bergstedt et al., 2012; Mattes, 2004; Moody et al., 2011). Therefore, ecotypes may be adapted to different growth temperatures, resulting in accelerated growth of the colder‐adapted siscowet, or slower growth for leans under experimental conditions (Levins, 1969; Stewart et al., 1983). Alternatively, growth of ecotypes in the wild was also probably influenced by differences in temperature and probably influenced growth observed by Hansen et al. (2016). Disentangling growth from temperature is often difficult, but variable growth rates are a common polymorphism across ecological gradients, and detectable in common garden experiments (De Villemereuil et al., 2016; Gardiner et al., 2010). Differences in growth can be a sign that ecotypes are on opposite sides of an ecological trade‐off (Grenier & Tallman, 2021). In this case, growth differences may be the result of differential investment into lipid stores between lean and siscowet. It is also important to recognize that our growth analysis was on fish only 8 years old which only provides a partial view of the total growth trajectory during the life span of these ecotypes (40+ years). In contrast, most von Bertalanffy analyses on wild populations, including Hansen et al. (2016), profile the growth through most of the life span in the population. Therefore, it is possible that the differences in growth parameters observed in our common garden is related to incomplete profiles of the growth history of these fish. Nonetheless, what can be concluded is that like lipid content, growth is fundamentally different between ecotypes.

Intermediate hybrid phenotypes

Environmental changes that increase gene flow can disrupt selective pressures isolating ecotypes and begin to homogenize a population through hybridization (Cenzer, 2016). Presently, conservation and management biologist are concerned that renewed gene flow among lake charr ecotypes in Lake Superior will result in the loss of distinct ecotype diversity (Baillie, Muir, Scribner, et al., 2016; Bronte et al., 2003; Ribeiro & Caticha, 2009). Our findings that hybrid crosses express intermediate lipid content suggests that siscowet and lean lake charr phenotypes are quantitative and additive in nature, and therefore could be lost through introgression. Siscowet and lean lake charr are generally believed to be reproductively isolated and spawn on different shoals, but given overlap in resident habitats, hybrid zones probably exist (Moore & Bronte, 2001). We cannot draw conclusions about the relative fitness of hybrids, but our study shows that hybrid crosses can reach maturity, may produce viable offspring, and express intermediate phenotypes of full crosses. Therefore, increased hybridization could facilitate the loss of distinct ecotypes and genetic diversity in Lake Superior.

Growth rate and size at age are often important factors that influence the fitness of offspring in fishes (Perez & Munch, 2010). Length and weight differences among all four crosses were much lower than differences in lipid content. However, while hybrid crosses expressed intermediate length and weight of full crosses in some years, they most often showed similar length and weight of the maternal lineage. This pattern suggests that maternal effects or inheritance may be contributing to length and weight, but not lipid content (Wolf & Wade, 2009). Maternal effects can have a substantial influence on Salmonid growth, especially early in life (March, 1991; Heath et al., 1999). Our results suggest that these effects may persist into the first several years of growth. Because full crosses did differ consistently in length and weight, the maternal effects we observed may be associated with differences in egg resource allocation between lake charr ecotypes. Our results indicate that maternal lineage may have an important influence on offspring fitness differences between lake charr crosses.

Polygenic regulation of trait variation

In many taxa, lipid content and growth differences are highly complex and likely to be regulated by many genes (Mak et al., 2006; Pomp, 1997). Our results support that this is also the case for lake charr ecotypes and indicate that lipid content and growth differences are most probably polygenically regulated. When the association of SNPs with lipid content and condition were compared across all crosses, using a standard GWAS approach, no chromosomes showed evidence of high‐effect loci or enrichment for trait‐associated markers. However, high family and geographic structure substantially lowered our power to detect associations using this approach (Santure & Garant, 2018; Sesia et al., 2021). By conducting RDA within‐cross design and requiring candidates to be identified independently in all four crosses, our ability to detect between‐cross differences was limited, but our probability of type II error was reduced. Using the markers identified with our within‐cross approach, polygenic PEA scores explained a significant amount of trait variance within each cross but not variance between crosses. These results indicate that many of the same loci which differentiate lake charr ecotypes may be involved with lipid metabolism in both ecotypes. Our results are supported by previous transcriptomic work showing that multiple genes in the liver involved in lipid metabolism and transport were expressed differentially between the ecotypes (Goetz et al., 2010, 2016). Finally, so far studies of wild populations have been unable to identify large effect loci differentiating ecotypes, and instead have hypothesized a polygenic mechanism is involved in phenotypic differences observed among ecotypes (Baillie, Muir, Hansen, et al., 2016; Perreault‐Payette et al., 2017). We analysed only a small subset of polymorphic sites across the lake charr genome and therefore our census of polygenic trait loci is incomplete, and it is possible that we missed some loci of large effect that explain variation among and within lake charr ecotypes (Barton et al., 2017; Szarmach et al., 2021). Furthermore, the unexplained variance between crosses may be controlled by undetected loci or epistatic relationships among genes (Lowry et al., 2017; Nosil et al., 2020). Additional analysis of these experimental groups in combination with wild fish could help identify loci differentiating in allele frequencies between crosses that may help to explain the variance in phenotype missed by our study. Our findings, however, do begin to describe the genetic contribution to key phenotypic traits that distinguish siscowet and lean lake charr ecotypes.

Gene ontology and QTL

Many of the candidate SNPs we identified were in or near genes suggesting that they may be linked to processes associated with lipid metabolism and growth in lake charr. Furthermore, the 95% CI around a QTL identified on Chromosome 38 contained at least three genes that coded for enzymes involved in the production or metabolism of fatty acids that could play a role in lipid metabolism, and we identified genes that are extremely important in the control of somatic growth and myogenesis including insulin growth factor 1 (Clemmons, 2012), and myogenic factors 5 and 6 (Asfour et al., 2018). This may suggest that the QTL we located is linked to other growth and energy metabolism traits that covary with lipid content. In prior studies with the parents of the F2s, liver transcripts involved in lipid binding and metabolism were differentially expressed between the ecotypes (Goetz et al., 2010, 2016), but the QTL region and candidate SNPs did not contain markers in any of those specific genes.

Conclusions and implications

Following the phenotypic variation of deep and shallow water lake charr ecotypes in a common garden experiment has helped to understand the ultimate mechanism maintaining phenotypic differences observed among lake charr in the wild. Long‐term experimental investigations, like ours, are essential to understand the dual role of environmental stochasticity and genetics in determining phenotypic variation in the wild. Through our investigation, we found that lipid content differences between lean and siscowet ecotypes appear to be heritable, suggesting that the mechanism for one of the key traits differentiating lake charr ecotypes in Lake Superior is primarily genetic and probably regulated by many genes that have an additive effect on phenotype. Meanwhile, standard metrics of growth, length, and weight appeared to be more likely maternally inherited, and ecotype specific. Our data supports previous studies that hypothesize that depth is an important ecological gradient in aquatic systems leading to resource polymorphisms and the evolution of sympatric ecotypes (Friedman et al., 2020; Lowry, 2012; Martin & Phennig, 2010). Many species experience similar environmental gradients and thus may also exhibit heritable differences across small spatial scales.

CONFLICT OF INTEREST

The authors declare no conflicts of interest for this study.

AUTHOR CONTRIBUTIONS

Peter T. Euclide drafted the manuscript and analysed the data. Andy Jasonowicz assisted preparing and processing molecular data and conducted laboratory measurements of phenotypic traits as well as assist with analysis. SS oversaw many components of the project, provided funding support, and assisted with drafting the manuscript. GF provided support housing and rearing of lake charr in the common garden and assisted in securing funding support. FG oversaw all components of the project, provided funding support, assisted drafting the manuscript, and initially conceptualized the project.

Supplementary Material

Click here for additional data file.