Temporal instability of lake charr phenotypes: Synchronicity of growth rates and morphology linked to environmental variables?

Abstract Pathways through which phenotypic variation among individuals arise can be complex. One assumption often made in relation to intraspecific diversity is that the stability or predictability of the environment will interact with expression of the underlying phenotypic variation. To address biological complexity below the species level, we investigated variability across years in morphology and annual growth increments between and within two sympatric lake charr Salvelinus namaycush ecotypes in Rush Lake, USA. A rapid phenotypic shift in body and head shape was found within a decade. The magnitude and direction of the observed phenotypic change were consistent in both ecotypes, which suggests similar pathways caused the variation over time. Over the same time period, annual growth increments declined for both lake charr ecotypes and corresponded with a consistent phenotypic shift of each ecotype. Despite ecotype‐specific annual growth changes in response to winter conditions, the observed annual growth shift for both ecotypes was linked, to some degree, with variation in the environment. Particularly, a declining trend in regional cloud cover was associated with an increase of early‐stage (ages 1–3) annual growth for lake charr of Rush Lake. Underlying mechanisms causing changes in growth rates and constrained morphological modulation are not fully understood. An improved knowledge of the biology hidden within the expression of phenotypic variation promises to clarify our understanding of temporal morphological diversity and instability.


| INTRODUC TI ON
A rapidly changing climate can have wide-ranging effects on organisms across ecosystems, which fosters a need to understand how ecosystems will respond to this variation in terms of structure and function (Montoya José & Raffaelli, 2010;Pacifici et al., 2015).
Contemporary climate change, including rapid increases in global temperatures, represents one of the most serious and current challenges to ecosystems, not only by threatening ecosystems directly (Norberg, Urban, Vellend, Klausmeier, & Loeuille, 2012), but also by contributing to cumulative and additive effects with other perturbations (e.g. industrial development, pollution, overhavest and nonnative species; CAFF, 2013;Poesch, Chavarie, Chu, Pandit, & Tonn, 2016). Ecosystems are mosaics of different habitats; climate change combined with abiotic and biotic variation across these habitats may lead to major eco-evolutionary responses (Grimm et al., 2013;Ware et al., 2019). Rapid biological responses to variation associated with environmental variation have already been detected at all levels, from individuals to species, communities and ecosystems (Heino, Virkkala, & Toivonen, 2009).
The importance of phenotypic variability has been emphasised in evolutionary and ecological population dynamics (Kinnison & Hairston, 2007;Schoener, 2011) because variation fuels evolutionary change (Stearns, 1989). Pathways through which phenotypic variation arises among individuals can be complex. Phenotypic variation can affect population dynamics such as through reproductive and mortality pathways (Bolnick et al., 2011). Furthermore, the magnitude of plasticity in the variation of trait expression differs among populations and ecotypes within a population (Skúlason et al., 2019). A conceptual framework to predict evolutionary and ecological consequences of climate change is currently limited by the scarcity of empirical data demonstrating phenotypic changes over time among individuals within ecosystems. The causes, patterns, and consequences of ecological and evolutionary responses to environmental variability need to be quantified across species, space and time.
Given rapid environmental changes are occurring within aquatic ecosystems, phenotypic differences among individuals must not be taken as negligible and ecologically inconsequential. Phenotypic variation is driven by switches along developmental pathways (West-Eberhard, 2003) that in some cases, can adjust immediately to variable environmental conditions (Japyassú & Malange, 2014).
Environmental factors can induce a component of variation that introduces fined-grained variation around coarse-scale temporal trends, resulting in year-to-year variation in phenotypes but not in genotypes-because genetic changes are not expected to be so responsive (Merilä & Hendry, 2014). Thus, phenotypic variability, both in means and variance, can provide an evolutionary scope for a population in the face of changing selection regimes by affecting population dynamics and probabilities of extinction (Chevin, Lande, & Mace, 2010;Johnson, Grorud-Colvert, Sponaugle, & Semmens, 2014;Reed, Waples, Schindler, Hard, & Kinnison, 2010).
Although models have often assumed stable ecological evolutionary equilibrium (Maynard, Serván, Capitán, & Allesina, 2019;Skúlason et al., 2019;Svanbäck, Pineda-Krch, & Doebeli, 2009), temporal variation in phenotypes among cohorts within an ecotype may not be ecologically trivial in a rapidly changing world. Lake charr are known to display intraspecific variation, mainly diversifying along a depth gradient, with shallow-versus deep-water ecotypes exploiting different prey resources within a lacustrine system (Chavarie et al., 2021). By selecting a case of co-existing shallow-and deep-water ecotypes of lake charr in Rush Lake, located at the southern edge of the species range ( Figure A1; , we examined the response of lake charr to environmental variation below the species level. Lake charr ecotypes in Rush Lake are likely to integrate diverse signals of climate variation across multiple trophic levels (Black, Matta, Helser, & Wilderbuer, 2013) because they exploit different prey items , invertebrates (deep-water ecotype) vs. foraging fish (shallow-water ecotype). Additionally, the effect of environmental variation, such as in temperature, on growth rates of aquatic organisms varies, especially across depths, in part because deep-water habitats are usually more environmentally stable (Jeppesen et al., 2014;Murdoch & Power, 2013;Thresher, Koslow, Morison, & Smith, 2007).
Considering that growth chronologies provide robust data sets for assessment of temporal and spatial environmental variation and its ecological consequences, variation across cohorts, ecotypes, species, and communities can be assessed (Black, 2009).
In this study, we tested whether phenotypic expression of two Rush Lake lake charr ecotypes remained stable or changed over time. We predicted that if a rapid change in phenotypic expression occurred and this change was associated with environmental variation, the shallow-water ecotype would display a greater magnitude of variation than the deep-water ecotype due to the shallow-water habitat being more responsive and sensitive to contemporary climate change. We measured phenotypic expression in terms of morphology and growth chronology. The objectives of our study were to (1) quantify the morphological variation between and within lake charr ecotypes over a ten-year period, (2) determine whether patterns of growth chronologies expressed by ecotypes were temporally synchronous with each other and associated with the phenotypic variation displayed within each ecotype; and (3) examine whether annual growth rates of lake charr ecotypes were related to environmental variation by using tree-ring cross-dating techniques.

| Assignment of lake charr ecotypes
Rush Lake is a small lake (1.31 km 2 ) that contains deep-water (>80 m) and is <2 km from Lake Superior .
Rush Lake provided the first documented example of sympatric lake charr ecotypes in a small lake, diverging along a depth gradient.
Two co-existing ecotypes of lake charr were found: a large streamlined-bodied shallow-water lake charr (lean ecotype) and a small plump-bodied deep-water ecotype (huronicus ecotype) ( Figure A1; Hubbs, 1929). Lake charr caught in 2007 were previously assigned to ecotype  and assignments for lake charr caught in 2018 used the same methodology (see Appendix B, Muir et al., 2014). Twenty sliding semi-landmarks and six homologous landmarks were digitized from each image to characterize head shape and 16 homologous and four sliding semi-landmarks were digitized from whole-body images to characterize body shape.
Morphological methods and programs are described by Zelditch, Swiderski, and Sheets (2012) and specific procedures were described by Chavarie, Howland, and Tonn (2013). All morphological measurements were size standardized by using centroid sizes (Zelditch et al., 2012).

| Temporal morphological variation between ecotypes and years
All data from 2007 and 2018 were combined to align all samples in the same shape space and partial warps for temporal morphological analyses between ecotypes and years. Principal component analysis (PCA) of body-and head-shape data was used to visualize morphological variation between and within lake charr ecotypes and years using PCAGen8 (IMP software). Canonical variate analyses (CVA) and validation procedures on body and head shape data were used to assess temporal differences within and between ecotypes using CVAGen (IMP software). Jacknife validation procedures included a test of assignment, with 1000 jackknife sets using 20% of our data as unknowns (Zelditch et al., 2012). Single-factor permutation multivariate analysis of variance (MANOVA) with 10,000 permutations of CVAGen was used to test whether body and head shape differed between and within (i.e. years) ecotypes. If MANOVA indicated differences among groups (α < 0.05), procrustes distance means were calculated for pairwise comparisons using TWOGROUP from the IMP software as post hoc tests (García-Rodríguez, García-Gasca, Cruz-Agüero, & Cota-Gómez, 2011). Procrustes distance for each pairwise comparison described the magnitude of difference between ecotypes and years. A bootstrapped Goodall's F statistic (N = 4900 bootstraps; full Procrustes based) was used to determine whether pairwise comparison differed. Allometric trajectories in body and head shape were compared between ecotypes and years by regressing PC1 scores (size standardized data) against centroid size (i.e., measure of size) (e.g., variation in developmental pathways can result in allometric trajectory patterns that can be parallel, divergent, convergent, or common; Figure A2; ; an allometric relationship occurred if the slope differed from 0. Relative body condition was compared between ecotypes and years in a 2-way analysis of variance (ANOVA), with main effects for ecotype and year and the interaction between ecotypes and years (Zar, 1999). To correct for size-related trends in body condition, relative body condition was defined as residuals from the power relationship between log 10 (W) and log 10 (TL) . If the ecotype*year interaction was significant, years were compared within ecotypes and ecotypes were compared within years in 1-way ANOVAs. To visualize the results, least-squares means (±SE) from the ANOVA were back-transformed from logarithms into original units of measure.

| Back-calculated length at age from otoliths: growth patterns displayed by ecotypes through time
Sagittal otoliths were used to estimate lake charr age and growth increments for fish sampled in 2007 and 2018. Otolith thin sections have been validated for age estimation of lake charr to an age of at least 50 years (for otolith preparation, see Appendix B; Campana, Casselman, & Jones, 2008). Otolith growth measurements can be used for several different purposes to gain ecological insight, but often need different analytical techniques to answer different questions. Towards this end, we used three analytical techniques and have provided a summary of the advantages and disadvantages of each (Table A1).
To determine whether patterns of growth chronologies expressed by ecotypes were temporally synchronous with each other, growth in length with age was modelled using a parameterization of the Von Bertalanffy length-age model (Gallucci & Quinn, 1979;Quinn & Deriso, 1999): The length-age model describes back-calculated length L t (mm) at age t (years) as a function of age at length = 0 (t 0 = years; incubation time of embryos from fertilization to hatching), early annual growth rate (ω = L∞ × K = mm/year; Gallucci & Quinn, 1979), theoretical maximum length (L ∞ = mm) and residual error (ε). Parameters ω and L ∞ were estimated using nonlinear mixed-effect models, package 'nlme' (R Core Team, 2016), with a fixed population effect (the average growth curve for the population from which individual fish were sampled), individual as random effect (growth curves for individual fish sampled from the population), and sex (male or female), ecotype (lean or huronicus) and sampled year (2007 or 2018) as fixed factors (to compare average growth curves between sexes, ecotypes and years; Vigliola & Meekan, 2009). Mixed-effects models are appropriate for modelling the within-group correlation of longitudinal, auto-correlated and unbalanced data, such as back-calculated growth histories (Pinheiro & Bates, 2000). Eight models of varying complexity were compared using AIC statistics (Burnham & Anderson, 2004): (1) ecotypes, sexes and sample years all included; (2) ecotypes and sample years both included; (3) ecotypes and sexes (4) sexes and sample years both included; (5) ecotypes only included; (6) years only included; (7) sexes only included; and (8) neither ecotypes, sexes nor years included.
Annual growth increments were modelled using a linear mixed-effects model (Weisberg, Spangler, & Richmond, 2010), wherein annual growth increments were modelled as a function of a fixed age effect (age of the fish when the increment formed), a random year effect (year in which the increment formed), a random fish effect (unique identifier for individual fish) and residual variation. The fixed age effect accounts for the fact that growth increments decline with age approximating a negative exponential curve. The random year effect reflects average increment width associated with each year of growth, after accounting for age effects (i.e., growth increment declines with age). The random year effect accounts for year-to-year environmental effects as random draws from a normal distribution, with a different draw for each year. The random fish effect accounts for fish-to-fish variation in growth as random draws from a normally distributed population with a different draw for each fish. The last source of variation is unmodelled residual variation. Differences between sexes (male or female), ecotypes (lean or huronicus), and sample years (2007 or 2018) were tested by including sex, ecotype and sample year as fixed effects.

| Otolith increment cross-dating: growth of lake charr in relation to environmental variation
In recent decades, advancements in dendrochronological techniques have been increasingly applied to sagittal otoliths, leading to novel insights on how broad-scale climate variation can impact both freshwater and marine systems (Black, Boehlert, & Yoklavich, 2005;Black et al., 2013;Matta, Black, & Wilderbuer, 2010). The use of dendrochronological methods (i.e. cross-dating techniques) ensures that specific growth annuli are assigned to an exact year (Black, Biela, Zimmerman, & Brown, 2012;Black et al., 2005Black et al., , 2013. In turn, this process enhances connection to common environmental signals across fish, which have syncronously limited growth in certain years. To assign otolith growth increments to exact calendar years, transverse sections of sagittal otoliths were aligned by calendar year and cross-dated visually using the list method (Yamaguchi, 1991) to confirm years when particularly large or small otolith increments would be expected. Thereafter, visual cross-dating was statistically confirmed using COFECHA software (Holmes, 1983). In using COFECHA, otolith time series with series intercorrelation values (i.e., Rbar) lower than 0.20 with the initial master chronology were removed and placed in a separate group that included more than one third of all fish sampled. In a previous study of lake charr otolith variation, Rbar values ranged from 0.42 to 0.97 (Black et al., 2013), supporting the assumption that otolith width series with Rbar <0.20 included anomalous inter-annual growth variation that did not match the initial master chronology. Otolith increment time series with Rbar <0.2 underwent a second round of cross-dating with COFECHA, separately from those that matched better with the initial master chronology. Because otolith increment data grouped together by ecotype without a priori knowledge of ecotype, all subsequent chronologies and analyses were conducted separately based on ecotype assignment from the morphological analyses described above.
Dendrochronological detrending methods generally attempt to remove growth variation and emphasize inter-annual variation in growth controlled by climate (Fritts, 2012). The first approach to detrending used the ARSTAN program (Cook & Krusic, 2014) to fit cubic splines of various rigidity based on fish age. Thereafter, autoregressive modelling was used to enhance inter-annual growth variability and the resulting indexed time series averaged within calendar years using a bi-weight robust mean. Thereafter, we undertook a second regional chronology standardization (RCS) approach known to better enhance low-frequency signals compared to the cubic spline method (Table A1; Briffa & Melvin, 2011). This detrending method divided each raw growth increment value by that expected from the mean growth increment for each ecotype and age combined ( Figure A3). These ratios were then multiplied by 100 and averaged within calendar years to yield a percentage change in growth for each year. In subsequent analyses of environmental influence on growth chronologies, data from each lake charr were truncated to feature only growth during young ages (age 1-3); which were excluded in earlier age-effect analyses but are a critical stage to phenotypic variation linked to environmental differences (Angilletta, Steury, & Sears, 2004;Georga & Koumoundouros, 2010;Ramler, Mitteroecker, Shama, Wegner, & Ahnelt, 2014). The approach employed for these comparisons corrected for age directly rather than using ARSTAN detrending (see Methods above and Table A1). We also limited the analysis to calendar years where each combination of ecotype and collection period included otolith data from at least seven fish. This constrained the calendar years investigated to 1986 to 2012 for huronicus ecotype and 1988 to 2010 for the lean ecotype.
Otolith increments, detrended with the ARSTAN program, were calculated as means within a year for both ecotypes and were initially compared against monthly resolution climate data for the corresponding year. Based on a priori knowledge of fish biology and lake-effect climate phenomena, temperature, precipitation and cloud cover were selected as environmental variables (Chavarie, Reist, Guzzo, Harwood, & Power, 2018;Voelker et al., 2019). Interpolated climate data (e.g. air temperature and precipitation) were obtained from ClimateNA version 5.6 software (Wang, Hamann, Spittlehouse, & Carroll, 2016). Cloud cover climate data from airports within 7 km of Lake Superior were obtained at daily resolution from the NOAA Great Lakes Environmental Research Laboratory (https://www.glerl. noaa.gov/) and summarized by month.
Pearson correlation values were calculated between each ecotype-specific growth chronology and monthly climate variables for the corresponding and previous two year (i.e., to detect lag effect).
After initial inspection of correlations between annual growth increments and monthly climate variables, the number of potential explanatory variables was consolidated into seasonal means, whereas winter to spring was defined as December through the next April for a corresponding year, summer was defined as June to August, and fall was defined as September to November. Autocorrelation was expected to be present in otolith increment time series and resulting chronologies due to year-to-year lags in growth owing to fat storage and subsequent metabolic withdrawals, skip spawning effects, and climate and climatic-effects on water temperature.
Robust assessments and modelling of autocorrelation on short time series are statistically impossible. Thus, we quantified what proportional weighting of climate data among the corresponding and two previous years produced the largest gains in Pearson correlations between otolith growth increment and seasonal climate data from an individual year to provide a window into how climate signals are incorporated into fish growth.
The influence of climate on otolith growth increment was assessed using forward selection multiple regression models, package 'lm' (R Core Team, 2016). More specifically: weighted temperature, precipitation, and cloud cover for each of four seasons as well as winter precipitation as snow (i.e., 13 total climate variables) were introduced iteratively to identify, for each ecotype separately, the variables that explained the most variation in otolith growth, which variables were significant (α < 0.05), and which combination resulted in higher Akaike Information Criterion (AIC) values (Burnham & Anderson, 2004) and were retained in the models. No models with more than two variables increased AIC values.

| Temporal morphological variation between ecotypes and years
In total, 107 lake charr were sampled, including 39 huronicus and 20 leans in 2007 and 27 huronicus and 21 leans in 2018. For both ecotypes, lake charr caught in 2007 had deeper bodies than lake charr caught in 2018 (Figure 1). The first two principal components explained 48.9% of the variation in lake charr body shape from Rush F I G U R E 1 PCA of lake charr body and head shape with percentage representing the variation explained by that component in (a) and (c), respectively, and CVA of lake charr body and head shape with 95% confidence ellipses delineating groups in (b) and (c) Lake ( Figure 1a). Body shape differed between years within each lake charr ecotype (CVA, Axis 1 ƛ = 0.015, p < 0.01 and Axis 2 ƛ = 0.29, p < 0.01; Figure 1b). Jackknife classification of body shape had a 54.3% rate of correct year and ecotype assignment (i.e., ecotypes and years as different factors). Body shape means differed between ecotypes and years (Permutation MANOVA, F = 11.7, df = 3, p ≤ 0.01), and the magnitude of these differences was slightly larger between ecotypes than between years. Pairwise body shape comparisons differed between ecotypes for both 2007 and 2018 (F tests; p ≤ 0.05; Figure A4). For lean and huronicus lake charr sampled in 2007, the Goodall's F was 21.1 and distance between means was 0.030 ± 0.0028 (SE), whereas for the lean and huronicus sampled in 2018, the Goodall's F was 18.2 and distance between means was 0.032 ± 0.0031 (SE). Pairwise body shape comparisons also differed between years for both ecotypes (F tests; p ≤ 0.05). For huronicus 2007 vs. 2018, the Goodall's F was 10.8 and distance be- For both ecotypes, lake charr from 2007 had deeper heads than lake charr from 2018 (Figures 1 and A4). The first two principal components explained 61.2% of the variation in lake charr head shape ( Figure 1c). Head shape differed between years within each lake charr ecotype (CVA, Axis 1 ƛ = 0.012, p < 0.01 and Axis 2 ƛ = 0.13, p < 0.01; Figure 1d). Jackknife classification on head shape had a 49.5% rate of correct year and ecotype assignment (i.e. ecotypes and years as factors). Head shape differed between ecotypes and years (Permutation MANOVA, F = 21.75 df = 3, p ≤ 0.01), and the magnitude of these differences was slightly larger between years than between ecotypes. Pairwise comparisons of head shape differed between ecotypes for both 2007 and 2018 (F tests; p ≤ 0.05, Figure A4). For lean and huronicus sampled in 2007, Goodall's F was 21.7 and distance between means was 0.040 ± 0.0035 (SE), whereas for the lean and huronicus sampled 2018, Goodall's F was 21.0 and distance between means was 0.048 ± 0.0044 (SE).
Pairwise head shape comparison also differed between years for F I G U R E 2 Allometric trajectories in body and head shape of lean and huronicus lake charr ecotypes from 2007 and 2018. PC1 scores of body shape (a) and head shape (b) are plotted against centroid size. Lean are represented by squares and huronicus by circles, whereas filled symbols are individual sampled in 2007 and empty symbols are lake charr caught in 2018. Only lean 2018 regressions were significantly different from 0, for both body and head shape (body shape: R 2 = 0.56, p ≤ 0.01; head shape: R 2 = 0.36, p ≤ 0.01) both ecotypes (F tests; p ≤ 0.05). For the huronicus, Goodall's F was 54.9 and distance between means was 0.063 ± 0.0036 (SE), whereas for the lean, Goodall's F was 38.6 and distance between means was 0.063 ± 0.0041 (SE). Allometric trajectories in head shape did not differ except for lean lake charr in 2018 (R 2 = 0.36, p < 0.01;

| Otoliths back-calculated growth: growth patterns displayed by ecotypes through time
Length at age of lake charr was best described by a single model that included both ecotypes and sample years (Table 1) were 44% more variable for lean than for huronicus ecotypes (i.e., growth varied more among years for leans than huronicus overall;

| Otolith increment cross-dating: growth of lake charr in relation to environmental variation
Inter-annual climate variation, particularly when including lagged effects, was correlated with fish growth for both ecotypes, as demonstrated by Pearson correlations regularly exceeding 0.2 ( Figure 6).
For both lake charr ecotypes and year-corresponding and lagged effects, annual otolith growth increments were positively correlated with summer air temperatures (except for the lean ecotype with corresponding year, slightly negative) and fall precipitation (i.e., more growth with warmer temperatures and more precipitation) and negatively correlated with summer and fall cloud cover (i.e., more growth with less cloud cover). In comparison, the direction of the relationship with winter to spring temperatures and precipitation as snow differed between the two ecotypes ( Figures 6   and A5). For each set of seasonal variables, inclusion of weighted climate data from previous years, to accommodate lagging effects, tended to strengthen correlations. Overall, based on the forward selection multiple regression models, the total amount of variation in otolith annual growth increments explained by climatic variables was greater for the lean ecotype than huronicus ecotype (R 2 = 0.56 vs 0.35; Table A3). These regression models confirmed that growth increments of huronicus ecotypes were most strongly associated with winter to spring temperatures and precipitation as snow and secondarily with summer temperature, whereas the lean ecotype was most strongly associated with winter to spring temperatures and previous fall precipitation (Table A3). Between ecotype differences in the relationship of annual growth increment with winter and spring temperatures and with precipitation as snow were then examined more closely, where regression analyses confirmed these differences (Figures 6 and A5). Specifically, differences in annual growth increments between ecotypes for any given year indicated F I G U R E 3 Relative body condition of huronicus and lean lake charr sampled from Rush Lake, in 2007 and 2018 that cold and snowy winters tended to favour growth for the huronicus ecotype whereas warmer winters with less snow favoured growth for the lean ecotype (Figures 6 and A5).

| D ISCUSS I ON
In this study, we demonstrated a rapid phenotypic shift that occurred for two lake charr ecotypes in a small lake. Within the last decade, a major decline in annual growth increments has occurred with a corresponding morphological shift in body and head shape for both lean and huronicus ecotypes. Even though the lean (shallow-water) ecotype displayed a greater annual growth variation over the years than the huronicus (deep-water) ecotype, both ecotypes displayed similar magnitude and direction of morphological changes in the most recent decade yielding an analogous 'sensitivity' in phenotypic change, meaning similar 'responsiveness' independent of habitats. We interpret these results to mean that the response threshold that determines individual sensitivity to a particular cue (i.e., environmental variables) caused similar phenotypic changes by individuals of both ecotypes (Baerwald et al., 2016;Moczek, 2010), in spite of differences in sensitivity between growth responses of the two ecotypes.
Our results concur with this latter mechanism because lake charr The effect of growth rate heterogeneity on morphological modulation (e.g., heterochonry, allometry; Klingenberg, 2014) has been observed to constrain or enhance morphological differences in several fish species (Heino, 2014;Jacobson, Grant, & Peres-Neto, 2015;Olsson et al., 2006). Exposure to different environmental conditions is known to have asymmetrical impacts on the stability of developmental pathways during early life stages (Lazić, Kaliontzopoulou, Carretero, & Crnobrnja-Isailović, 2013;Robinson & Wardrop, 2002

TA B L E 2
Tests of fixed effects for differences between lake charr morphs (huronicus or lean), sample years (2007 or 2018), and sexes (males or females) from a linear mixed-effects model of annual sagittal otolith growth increments as a function of a fixed age effect (age of increment formation), random year effects (year of increment formation) and random fish effects (fish-to-fish variation in growth) sampled in Rush Lake F I G U R E 5 Annual growth increments (mm; random year effects from a linear mixed-effects model that also included fixed age effects and random fish effects; Weisberg et al., 2010) by calendar year for huronicus and lean lake charr ecotypes sampled from Rush Lake, in 2007 and 2018 Allometry, the shape variation associated with size variation (Zelditch et al., 2012) is also seen as a canalized process and an interacting agent that can limit morphological variation (Klingenberg, 2010). At the very least, some of our results represent a case of plastic allometry for the lean ecotype in 2018. Considering that the lean body condition did not markedly differ between sampling years and was within the species range values (Hansen, Guy, Bronte, & Nate, 2021), we are confident that morphological differences were not due to starvation of individuals.
When compared to lake charr ecotypes from other lakes, the genetic diversity and divergence of lake charr in Rush Lake is low for both ecotypes . This low genetic diversity and divergence favour the hypothesis that phenotypic variation is the result of phenotypic plasticity rather than genetic adaptations (although rapid genetic change cannot be excluded). Epigenetically mediated biological complexity is known to be an important process to tailor phenotypic reaction norms (e.g., linear and nonlinear)  Svanbäck et al., 2009). Such spatial and temporal periodicity can occur at different scales (e.g., for temporal periodicity: seasonal, inter-annual, decadal) and be driven by fluctuations in abiotic (e.g., temperature, precipitation, light and nutrients) and biotic processes (e.g. growth, reproduction and trophic interactions), thereby shaping species behaviour and 'rewiring' food webs (Bartley et al., 2019;McMeans, McCann, Humphries, Rooney, & Fisk, 2015).
In this study, some aspects of body and head shape shifts may have been mediated by changes in growth rates and body condition, which likely, were induced by variation in food availability and abiotic conditions. Body condition can affect body shape in fish, and often reflects bulkiness of individuals, and consequently body depth (Borcherding & Magnhagen, 2008;Jacobson et al., 2015;Olsson et al., 2006;Svanbäck et al., 2017). Although both ecotypes were subjected to annual growth rate declines, the huronicus ecotype was more affected by body condition changes than the lean ecotype, which might suggest physiological differences between ecotypes in the dynamics of energy processing (e.g., metabolism and reproduction) and storage (Goetz et al., 2013). For example, in 2018, 30.8% of the huronicus females were in a resting reproductive stage compared to only 8.3% of lean females. In contrast to females, all males in 2018 of both ecotypes were reproductively mature (Unpublished data). Although no sexual dimorphism occur in lake charr (Esteve, McLennan, & Gunn, 2008), differences in energetic requirements associated with reproductive output appear related to the apparent skipped spawning patterns observed by ecotype and sex. Many examples exist of taxa where population growth fluctuated over time as a result of variation in resource levels, often were influenced by environmental changes (Ohlberger, 2013;Persson & De Roos, 2003). Given that fish otoliths offers a unique broad comparative tool to link abiotic factors as driving size changes (Gardner et al., 2011), we were able to detect similarities and differences of annual growth rates of lake charr ecotypes correlated to environmental variation. In our study, cloud cover was the main environmental variable that had steadily decreased over the same time period that lake charr growth declined (except for early stage) and morphology shifted in Rush Lake (Figure 7). Environmental heterogeneity is thought to have stronger effects on morphology at early life stages (Johnson et al., 2014;Morris, 2014;Ramler et al., 2014), suggesting that the effects of the cloud cover on lake charr might have been more significant at age 1-3 years than during later stages of life. The effect of cloud cover on growth at age 1-3 was stronger in the lean than the huronicus ecotype, which might explain why allometry was detected only for the lean ecotype in 2018. The relationship between lake charr annual otolith growth and cloud cover could be related to how solar irradiance can co-vary with other climate variables that may affect fish growth (Poesch et al., 2016;Reist et al., 2006).
Additionally, when higher temperatures are accompanied by suitable net addition of food ration (e.g. from direct and indirect effects of temperature and precipitation factors), increases in growth could be expected up to the point of the optimum metabolic temperature of the species (Elliott & Elliott, 2010;Elliott & Hurley, 2003).
Response to winter climate appeared to vary between ecotypes but it is unclear why the deep-water ecotype would show higher annual growth increments in years with warm winter temperatures and low snow cover (via direct or lagging effects) and the shallow-water F I G U R E 7 Otolith increment growth variation among ecotypes and collection periods for ages 1-3 plotted versus regional cloud cover from July to September for the same year without weighting or lag-effects included in (a) and (b). Data were detrended using regional chronology standardization (see Methods and Table A1). Within a given year, ecotype and collection period, otolith growth data for earlystage growth were averaged across ages 1-3. In (c), reductions in regional cloud cover by month from 1983 to 2017. Shown here are the predicted values for 1983 and 2017 using linear regressions fitted to inter-annual cloud cover data by month across the same time period. These two years represent endpoints for which otolith data were replicated enough to allow chronology construction. Cloud cover data were from airports within 7 km of Lake Superior since the lake itself strongly controls nearby weather (see Voelker et al., 2019) ecotype would have higher annual growth increments in years with cold winter temperature and high snow level. These results could be explained in part by reduced habitat partitioning between ecotypes during winter, along with an increase of intraspecific competition, affecting energy storage (Amundsen, Knudsen, & Klemetsen, 2008).
Another explanation, which is not exclusive of the previous one, could be that each ecotype's prey types, density, and quality (e.g., time response to environmental variable) are modified differently by lagging effects from winter environmental conditions (e.g., ice cover duration, ice and snow thickness). The lean ecotype is known to feed on forage fish whereas the huronicus ecotype mainly feeds on the invertebrate Mysis . Invertebrates can differ in response time and magnitude to environmental changes compared to forage fishes (Heino et al., 2009, Wrona et al., 2016. This hypothesis of a differential response of prey items to environmental conditions seems to be strengthened by correlations when weighted climate data from previous years were included with lake charr growth. This result might be expected for an organism in which growth rates may draw on a mixture of recently acquired and stored resources or where climate variables in one year may affect the abundance and composition of prey in subsequent years. For the most part, winter environmental conditions can play an essential role in ecological and evolutionary processes that define life-history characteristics (e.g., somatic growth, size and age at maturity, reproduction investment and longevity) of lacustrine species (Shuter, Finstad, Helland, Zweimüller, & Hölker, 2012).

| CON CLUS ION
One assumption often made in relation to intraspecific diversity, mostly tested experimentally or modelled, is that a stable or predictable environment interacts with underlying variation in expression of phenotypes (Skúlason et al., 2019;Wagner & Schwenk, 2000). In our study, the magnitude and direction of the observed phenotypic shift in both annual growth and morphology over a single decade were consistent for each ecotype and suggested similar pathways to which phenotypic variation was expressed. The degree of phenotypic variation that occurs within an ecotype theoretically depends on the relative strength and timing of mechanisms that drive phenotypic change (Wood et al., 2020). In our case, the observed phenotypic shift was relatively small, but nonetheless, detectable (i.e., cryptic eco-evolutionary dynamics; Kinnison, Hairston, & Hendry, 2015).
Several questions arise from our results, but one of interest is the organism's capacity for phenotype acclimation to a changing environment (via phenotypic plasticity and adaptation; Gorsuch, Pandey, & Atkin, 2010;Huey & Berrigan, 1996). Was this phenotypic shift an isolated event or does this type of change occur frequently in this lake and elsewhere? Answers to the question of phenotypic acclimation and the frequency of its occurrence within and among systems would require long data sets collected over multiple decades and would help to fill an important knowledge gap about nonequilibrium population dynamics affecting evolutionary dynamics.
Mechanisms that connect annual growth increments with morphological modulation are not fully understood (Olsson et al., 2006(Olsson et al., , 2007Svanbäck et al., 2009); however, the biology underlying phenotypic variation can have major implications for populations responding to climate change. Multidimensional phenotypic variability and its influence on patterns of population dynamics is a relatively poorly studied phenomena (Westneat et al., 2015), but individual and population resistance and resilience to climatic changes may depend on this variability (Johnson et al., 2014). The similarity in phenotypic response expressed by both ecotypes raises the question whether organisms in small lakes are more vulnerable to climate change than those in large lakes. Small lakes generally sustain a higher degree of habitat coupling (e.g., littoral-pelagic; Dolson, McCann, Rooney, & Ridgway, 2009;Schindler & Scheuerell, 2002), which is critical to food-web dynamics. Thus, the degree of habitat coupling found in each freshwater ecosystem might translate to its degree of vulnerability to climate change. Field studies, such as ours, that focus on temporal phenotypic instability within an aquatic ecosystem promise to clarify our understanding of how the interplay among phenotypes, trophic dynamics, and environmental context influences both ecosystem and evolutionary processes (Ware et al., 2019).

ACK N OWLED G EM ENTS
We thank the Huron Mountain Club for access to their lands, hous-

CO N FLI C T O F I NTE R E S T
The authors declare that they have no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data sets supporting the conclusions of this article are included within the article. Raw data will be available on Dryad (https://doi.  Figure 5, Figure 6, Table A2 Cross-dated and ARSTAN-detrended Conventional approach for identifying climate drivers of inter-annual variation in long-lived organisms such as trees Short age span of fish excludes identification of decadal or greater environmental drivers (i.e. 'segment length curse'). Some loss of data when otolith time series do not cross-date Figure 7, Table A3, Figure A5.

Cross-dated and RCS-detrended
Alternative approach for identifying climate drivers that preserves decadal and longer signals Not as efficient at accurately extracting inter-annual variation compared to ARSTAN-detrending Figure 8, Figure A3 TA B L E A 2 Difference in annual growth increments (mm; random year effects from a linear mixed-effects model that also included fixed age effects and random fish effects; Weisberg et al., 2010) by calendar year for huronicus and lean lake charr ecotypes sampled from Rush Lake, in 2007 and 2018 (see Figure 6)

F I G U R E A 1
Lean-like (a) and huronicus (b) lake charr Salvelinus namaycush ecotypes of Rush Lake. Illustration by P. Vecsei, after . The lean ecotype has a long head, long maxillae and a posterior eye position, which are all characteristics for piscivorous feeding (Janhunen, Peuhkuri, & Piironen, 2009;Keeley, Parkinson, & Taylor, 2007;Proulx & Magnan, 2004). The huronicus ecotype, with smaller gape and higher eye position than the lean ecotype, appears adapted for lowlight vison and as a vertical migrating predator feeding on the invertebrate opossum shrimp (Mysis spp.) as its main prey (Hrabik, Jensen, Martell, Walters, & Kitchell, 2006;Muir et al., 2014) F I G U R E A 2 Possible pattern of allometric trajectories between two putative ecotypes, figure is modified from . Allometric trajectory patterns can be (a) parallel, (b) divergent, (c) convergent or (d) common and minimum and maximum water depth were recorded for each net set. All fish caught were photographed in lateral view (Muir et al., 2014), and sagittal otoliths were removed for analysis of age and growth.

Assignment of lake charr ecotypes
Ecotypes were assigned to each individual using a combination of Bayesian cluster analyses using head and body shape information Disagreements between model and visual assignments were settled using decision rules described by Muir et al. (2014).

Age assignment and growth increments from otoliths
One otolith from each fish was embedded in epoxy. A Buehler Isomet 1000 Precision Saw was used to remove a thin transverse section (400 µm) containing the nucleus perpendicular to the sulcus.
Sections were mounted on glass slides and polished. Digital images of otolith sections were captured for age and growth assessment.
Criteria for annulus demarcation followed those of Casselman and Gunn (1992). Age estimates were used to inform demarcation of growth increments, measured from the nucleus to the maximum ventral radius of the otolith, and radial measurements at each annulus were used to back-calculate length at age using the biological intercept back-calculation model (Campana, 1990). The biological intercept (sagittal otolith radius = 0.137 mm; age-0 lake trout length = 21.7 mm; Hansen et al., 2012) was based on equations describing