Testing the predictability of morphological evolution in contrasting thermal environments

Gaining the ability to predict population responses to climate change is a pressing concern. Using a “natural experiment,” we show that testing for divergent evolution in wild populations from contrasting thermal environments provides a powerful approach, and likely an enhanced predictive power for responses to climate change. Specifically, we used a unique study system in Iceland, where freshwater populations of threespine sticklebacks ( Gasterosteus aculeatus ) are found in waters warmed by geothermal activity, adjacent to populations in ambient-temperature water. We focused on morphological traits across six pairs from warm and cold habitats. We found that fish from warm habitats tended to have a deeper mid-body, a subterminally orientated jaw, steeper craniofacial profile, and deeper caudal region relative to fish from cold habitats. Our common garden experiment showed that most of these differences were heritable. Population age did not appear to influence the magnitude or type of thermal divergence, but similar types of divergence between thermal habitats were more prevalent across allopatric than sympatric population pairs. These findings suggest that morphological divergence in response to thermal habitat, despite being relatively complex and multivariate, are predictable to a degree. Our data also suggest that the potential for migration of individuals between different thermal habitats may enhance nonparallel evolution and reduce our ability to predict responses to climate change.

Understanding whether populations evolve in a predictable manner when exposed to similar environments is crucial for understanding adaptation.Studies on a wide range of taxa (insects: Nosil et al., 2002;fishes: Bernatchez et al., 2010;birds: Mundy, 2005; and mammals: Hoekstra, 2006) have shown that different populations in similar environments tend to evolve similar phenotypes (Losos, 2011).Such consistent patterns are unlikely to be attributable to random processes and can arise due to natural selection, developmental bias, or a combination of the two (see Parsons et al., 2016;Uller et al., 2018).Thus, evolutionary forces can often drive population divergence in a way that is predictable and repeatable (Brakefield, 2006;Oke et al., 2017;Schluter & Nagel, 1995.Identifying the occurrence of shared patterns of evolution is often indicative of adaptation.Recently, this has been shown to occur in response to anthropogenic stressors and could be especially valuable for providing insights in light of global environmental change (de Amorim et al., 2017).This is because the ability to predict general population responses to anthropogenic change is of central importance for the planning of management and conservation efforts.In the coming decades, climate change will arguably pose the most significant threat to biodiversity.Rising temperatures are already altering abiotic and biotic environmental conditions and imposing novel selection pressures as well as developmental conditions (Campbell et al., 2017;Crozier & Hutchings, 2014).Ectotherms, such as fishes and reptiles, are particularly vulnerable because of their high sensitivity to temperature changes (Zuo et al., 2011).Consequently, there is now a pressing need to understand the scope for populations to respond to climate change.
Climate change will occur over decades, but the timescale of many studies on the impacts of increased temperatures focuses only on plastic (within-generation) responses.Studies directly investigating long-term evolutionary responses to climate change are still relatively scarce (Crozier & Hutchings, 2014).While contemporary evolution has been demonstrated in response to a number of environmental gradients, there are few available studies examining evolutionary responses to different temperature treatments over several generations beyond those conducted in laboratory settings (e.g., Alton et al., 2017;Barrett et al., 2010).Such experiments can usually only examine the direct effects of temperature, but under natural conditions, changes in temperature will be accompanied by changes in various ecological factors, such as food availability, parasitism, and predation pressure (Crozier & Hutchings, 2014).We, therefore, propose that examining natural populations inhabiting contrasting thermal environments provides a more powerful approach to understanding and predicting population responses to increasing temperatures.
To this end, we took advantage of a "natural experiment" in Iceland, where freshwater populations of threespine sticklebacks (Gasterosteus aculeatus) are found in waters warmed by geothermal activity (warm habitats), adjacent to populations in ambient-temperature water (cold habitats).This unique study system is relatively new to science (Jutfelt, 2020) and provides repeated and independent examples of populations inhabiting long-term contrasting thermal environments over a small geographic scale, thereby avoiding the confounding factors associated with latitudinal or elevational comparisons.Although there are some differences in water chemistry between thermal habitats, most variables are similar, suggesting that temperature could play a key role in potential divergence (see online supplementary material, Table 1).Another notable attribute of these study systems is that while most of these warm and cold habitats are in separate water bodies with geographic barriers (allopatric), some are found in different parts of the same water body with no physical or geographic barriers to migration (sympatric).Movement of individuals, and thus gene flow, is more possible between sympatric populations relative to allopatric populations, while in allopatry movement between habitats would be extremely rare (e.g., the dropping of fish by bird predators).This allows us to examine whether connectivity might influence the magnitude and/or direction of divergence between warm and cold habitats.Lastly, the age of these warm habitats, and hence the maximum time these populations have experienced elevated temperatures, ranges from decades to thousands of years (Table 1).These different timescales should make it possible to examine whether populations exposed to a warm environment for a relatively short contemporary timescale have diverged to the same extent as much older populations.
Here, we focus on providing a formal examination of thermally driven divergence and evolution in morphology.Specifically, establishing whether morphology has diverged between thermal habitats, and in turn, whether this has a heritable basis will provide fundamental knowledge about the impact of temperature on natural populations, while also establishing trends for a broader investigation.Morphology provides an inroad toward this, as it frequently contributes to fitness by influencing reproduction, foraging ability, and swimming performance (Rowiński et al., 2015).Morphological variation is also related to sexual selection and reproductive isolation and can therefore contribute to population divergence and speciation (Head et al., 2013).Previous research has shown that morphology in fishes often exhibits similar patterns of divergence across populations in response to common environmental conditions (Cooper et al., 2010;Jastrebski & Robinson, 2004;Stuart et al., 2017).While little is known about the effects of thermal conditions on body shape, an emerging trend is that warmer waters can be associated with increased body depth (Fruciano et al., 2011;Piñeros et al., 2015;Rowiński et al., 2015).More proximately, the temperature has been shown to influence the initiation of bone ossification, occurring at an earlier stage at warmer temperatures, with accompanying increases in bone remodeling that can have fitness consequences (Campbell et al., 2021).From laboratory studies, we also know that rearing temperature directly influences the development of body shape in fishes through plasticity (Campbell et al., 2021;Ramler et al., 2014;Sfakianakis et al., 2011).Yet, it is still unknown whether such morphological responses to temperature can result in the evolution of shape over generations, and whether such outcomes could be predictable across populations.Evidence for predictable evolutionary responses to conditions found in natural habitats with elevated temperatures could be extremely valuable, but these study systems will also allow us to test for a range of outcomes.For example, while population pairs may display shared patterns of divergence between thermal habitats, they may not necessarily show the same magnitude of divergence.Hence, we considered (a) whether divergence between thermal habitats was present across populations, (b) whether thermal divergence occurred across similar trajectories, and (c) whether the degree of thermal divergence differed across population pairs.Also, in correspondence with habitat connectivity, we tested (d) whether allopatric divergence differed from sympatric divergence.Lastly, to determine whether patterns of morphological variation observed in wild-caught sticklebacks were heritable, we conducted a common garden experiment by breeding fish from warm and cold habitats and rearing their offspring at a common temperature.By addressing these aspects in a natural system, we can demonstrate that evolutionary responses in morphology are likely to take place and are potentially predictable in the conditions provided by elevated temperatures.

Collecting sticklebacks
We used unbaited minnow traps to collect adult threespine sticklebacks from six warm-cold population pairs in Iceland in May and June of 2016 (Table 1, Figure 1).Three of these population pairs were allopatric (designated A1-3), meaning that the warm and cold habitats were in neighboring but separate water bodies (Table 1).The other three population pairs were sympatric (designated S1-3), meaning that the warm and cold habitats were found in the same water body with no physical barriers between them (Table 1).Indeed, population genomic analyses indicate extensive gene flow between thermal habitats in sympatric pairs (Costa et al., unpublished data).However, despite this, there is evidence of physiological divergence between thermal habitats in both sympatry and allopatry, with warm-habitat fish having a reduced standard metabolic rate (Pilakouta et al., 2020).
All cold habitats except A1 have existed since the last glacial period about 10,000 years ago (Einarsson et al., 2004).The A1 cold habitat appears to be a man-made trench likely dug during housing construction in the 1940s.There is also some variation in the age of the warm habitats (Table 1).The A1, S2, and S3 warm habitats originated 50-70 years ago and are fed by excess hot water run-off from nearby residences that use geothermal heating.The remaining warm habitats have been naturally heated by geothermal vents for over 1,000 years (Einarsson et al., 2004;Hight, 1965).
A subset of sticklebacks (n = 331) caught in minnow traps was immediately euthanized using an overdose of phenoxyethanol and preserved in 10% buffered formalin (Table 1).In addition, approximately 100 sticklebacks from each of the eight sampling locations (i.e., four warm-cold population pairs) were kept in temporary holding tanks at the Hólar University before being transported to the University of Glasgow to be used for breeding in our common garden experiment (see below).

Transport of study animals and animal husbandry
We fasted sticklebacks for 48 hr to minimize the build-up of ammonia in the transport water.On the day of shipping, we placed approximately 100 sticklebacks from each population in 100-L polyethylene bags containing 25 L of water.Air was removed from the bags and replaced with pure oxygen.Bags were sealed and placed inside insulated Styrofoam shipping boxes to minimize temperature fluctuations during transport.The fish were in transit for approximately 72 hr before arriving in Glasgow.No mortality was observed during transport.
Once the fish arrived at the University of Glasgow, they were kept at densities of 10-15 individuals per 10-L tank in a common recirculation system at 15 °C.This intermediate temperature is close to the maximum temperature experienced by fish in cold habitats in the summer and the minimum temperature experienced by fish in warm habitats in the winter (Pilakouta et al., 2020).All tanks contained plastic plants as shelter and air stones to oxygenate the water.Fish were fed ad libitum twice a day with a mixture of frozen bloodworms, Mysis shrimp, and Daphnia.They were kept at a 12-hr light:12-hr dark photoperiod.

Common garden experiment
We carried out a common garden experiment to determine whether morphological variation between fish from warm and cold habitats was heritable.For this experiment, we bred wild-caught sticklebacks from warm and cold habitats of two allopatric population pairs (A1 and A2) and two sympatric population pairs (S1 and S2).Gravid females and males displaying breeding colors were euthanized with an overdose of benzocaine and used for in-vitro fertilization (Barber & Arnott, 2000).We dissected males to remove their testes, which were then macerated to release the sperm.We stripped egg clutches from females by applying gentle pressure to the abdomen (Barber & Arnott, 2000).We then performed in-vitro fertilization in Petri dishes by mixing the sperm and eggs from one male and one female to create full-sib families.Fertilized embryos were placed in mesh baskets submerged in well-aerated 18 °C (±0.5 °C) water with methylene blue (2.5 µg/ml) until hatching.
These F1 generation stickleback larvae were fed with newly hatched HUFA-enriched Artemia salina nauplii, microworms, and powdered food (ZM100 and ZM200 fry food, ZM Systems, Twyford, UK) until large enough to eat pelleted food (Microstart, EWOS Ltd, Surrey, UK) at a standard length of approximately 2 cm.At this stage, fish were transferred to 10-L tanks and kept at standardized densities of 15-20 individuals to prevent large differences in growth rate between families.They were maintained at a constant water temperature of 18 °C (±0.5 °C) from the embryonic stage to adulthood.About 12 months after hatching, we euthanized 320 F1 individuals using an overdose of benzocaine and preserved them in 10% buffered formalin (n = 40 per sampling location from at least five full-sib families; Table 1).

Specimen preparation
All preserved specimens were bleached and cleared to remove skin pigmentation and make the body translucent (Potthoff, 1984).They were then stained with Alizarin Red S to emphasize bone morphology and were stored in 75% glycerol until the excess stain was removed.Individual specimens (n = 331 wild-caught sticklebacks, n = 320 F1 sticklebacks) were photographed on their left side with a Canon EOS 1100D digital camera (Canon Inc, Tokyo, Japan).All photographs included a scale and were taken from a fixed distance and angle using a copy stand.

Linear measurements
We used photographs of wild-caught and F1 sticklebacks to measure pectoral fin length and dorsal spine length, which are related to swimming and defense, respectively (Drucker et al., 2005;Hoogland et al., 1956).These traits were of particular interest as they can be directly attributed to bone length, and temperature can affect bone growth in fishes (Campbell et al., 2021).Using the software program tpsDig2 (Rohlf, 2018), we placed landmarks on the base and tip of the longest fin ray (pectoral fin length) and on the base and tip of the first and second spine (dorsal spine lengths).We then calculated raw interlandmark distances based on the Pythagorean theorem using the "Tradlength" routine within CoordGen8 (Sheets, 2018).Because they belonged to articulated structures, these landmarks were only used to obtain linear distances and were excluded from multivariate body shape analyses.After checking for homogeneity of slopes using an analysis of covariance, we regressed each linear distance against geometric centroid size and used the residuals in a separate analysis of variance (ANOVA) models (described below) for each linear trait.Geometric centroid size represented overall body size and was calculated as the square root of the sum of all squared distances to the geometric center of each specimen.

Quantifying body shape variation
To quantify morphological shape in wild-caught and F1 sticklebacks, we placed 22 anatomical landmarks on each image using tpsDig2 (Rohlf, 2018) to assess variation in the left lateral view using a geometric morphometric approach (Figure 2).
To minimize the effects of size and orientation across individuals, we then performed a single Generalized Procrustes Analysis on all specimens (wild-caught and lab-reared simultaneously) using the gpagen function within the "geomorph" package (Adams et al., 2020) using the R programming language (version 4.2.1, (R Core Team, 2018).This process superimposed landmark configurations to minimize the sum of squared distances between corresponding landmark configurations by scaling, rotating, and translating specimens in relation to their geometric center.
Following this, we minimized the potential effects of curvature.Due to the soft bodies of fish, preservation processes can introduce uninformative variation in body curvature despite efforts to minimize such issues during the photographing stage.This variation was visually identified through a Principal Components Analysis using gm.prcomp, which identified that the first principal component for each of the wild and F1 lab-reared fish was primarily influenced by body curvature.Therefore, multivariate regressions were performed using shape coordinates from each of the wild and F1 lab-reared fish on their corresponding PC1 scores using procD.lm.To produce "straightened" landmarks, the predicted forms with the least bent shape were determined using shape.predictor to which the residuals from the multivariate regressions were added.

Statistical analysis
Statistical analyses were run using R version 4.2.1 (R Core Team, 2018) and figures were generated using the "geomorph" (Adams et al., 2020) and "ggplot2" packages (Wickham, 2016) unless otherwise noted.

Linear measurements
We used ANOVA models for each linear trait to test the effects of population pair, thermal habitat, and their interaction on their residuals.The main effect of population pair summarized the properties unique to different replicates (Bolnick et al., 2018).The main effect of thermal habitat measured the extent to which thermal divergence is shared across sampling locations (Bolnick et al., 2018).The population pair × thermal habitat interaction indicated how the direction and magnitude of thermal divergence varied among population pairs, accounting for unique types of divergence.To determine the partial variance explained by each factor and interaction, we used the "heplots" R package to calculate partial eta squared (η 2 ) values (Fox et al., 2007).Lastly, Tukey's HSD posthoc tests were performed using TukeyHSD to identify which population pairs were responsible for potential differences in linear traits between thermal habitats.

Body shape variation
To test whether thermal habitat affected body shape we applied models based on residual randomization in permutation procedures (RRPP) using the "RRPP" package.This method is especially useful for high-dimensional data where the number of variables exceeds the number of observations (Collyer & Adams, 2019).We applied RRPP models using ordinary least squares and type 1 sums of squares and cross products.Our first RRPP model included all wild-caught and lab-reared fish and used population pair, thermal habitat (warm and cold), and allopatry/sympatry status (nested within a population) as factors.Also, to account for potential allometric effects we included geometric centroid size as a covariate within this model.To obtain the simplest and best supported RRPP model we compared two models, one including all factors against one excluding allopatry/sympatry status.Using a nested ANOVA within RRPP we tested their ability to explain shape variation.Both models were run with 1,000 permutations and ANOVA was performed on each using random distributions of F-statistics to calculate z-scores, r-square, and p-values.While these models could assess overall trends by including all fish, it was important to identify whether thermal divergence was present within each population pair.Therefore, to assess divergence within each pair we used RRPP models (one for each of the wild, and when available, lab-reared fish) including thermal habitat and geometric centroid size to explain shape variation.Similar to the first model, all population pair-specific models were run with 1,000 permutations with an ANOVA performed to provide summary statistics.
For each population pair, we then visualized differences between geothermal and ambient fish through the generation of deformation grids.This was carried out by reciprocally comparing the consensus configuration of each warm and cold group from each of our population pairs, including both wild-caught and F1 lab-reared fish.All deformation grids were magnified by a factor of 3 to accentuate differences and to aid in interpretations of shape, and were produced using the plotRefToTarget function within "geomorph."

Direction and magnitude of thermal divergence in body shape
While the main effect of thermal habitat in our best RRPP model (that included all wild and lab-reared fish) could indicate evidence of shared patterns of divergence, the proximate shape components (including the direction and magnitude of divergence) required additional steps to reveal.Specifically, we tested for evidence of parallel phenotypic divergence and examined the magnitude of divergence at a finer pairwise scale between the divergence trajectories of population pairs.This trajectory analysis was performed using a trajectory analysis within the "RRPP" package using fitted values from our best model that included all study fish (Adams & Collyer, 2009;Collyer & Adams, 2019).
From a strict perspective, "parallelism" in phenotypic divergence should be represented by completely parallel divergence trajectories (i.e., 0° difference in trajectory).However, in most biological systems undergoing divergence this level of stringency is likely unrealistic.Within our trajectory analysis, parallelism was the null hypothesis, meaning that significant differences in trajectory could be used to reject parallelism but a lack of difference could not be used to accept it.This test was also based on confidence limits meaning that angles considerably larger than 0° could be non-significant (i.e., suggesting parallelism).Therefore, we arbitrarily defined parallelism as the presence of divergence trajectories that differed no more than 15°, a strict criterion, but in line with what might be expected biologically.Nonetheless, beyond our stricter criterion, this test did provide a means to more directly assess how similar the direction of divergence was between population pairs.We used summary.trajectory.analysis to quantify multivariate trajectories defined by the factor levels within our RRPP model.This function compared trajectory attributes between each population pair to test for significant differences in the magnitude of divergence as well as divergence in trajectory orientations (i.e., the aforementioned tests of parallelism).Magnitudes of divergence were calculated as Procrustes distances between the least square mean shape of populations, while trajectory orientations were derived from a subtraction of mean coordinates with the product used to test for correlations between divergence vectors.Both the magnitude differences between trajectories and trajectory correlations were then calculated from each of the 1,000 residual randomization permutations following Collyer and Adams (2019).We performed these comparisons for all combinations of our population pairs (i.e., 15 pairwise comparisons for wild fish, four pairwise comparisons between wild and F1 fish from the same pair, and six pairwise comparisons among F1 fish pairs).
Finally, because our trajectory analysis used parallelism as the null we performed an additional test, where we instead tested for significant similarity in vector orientation.
Here, the null hypothesis was that vector orientations would be no more similar than expected by chance.To base these tests on equivalent data we extracted the fitted values of our best RRPP model which included all wild and lab-reared fish and calculated the least square coordinate means for geothermal and ambient fish within each population pair.Divergence vectors were then calculated as before, by subtracting coordinates to obtain a product TestOfAngle within the "GeometricMorphometricsMix" package (https://www.fruciano.org/software/).Briefly, the main purpose of this function was to test for similarity in trajectory orientation on the basis of a closed-form formula for the area of a hyperspherical cap (Li, 2011).The area of the cap divided by the area of the entire hypersphere indicated the probability that a random vector drawn from a uniform distribution form an angle with a fixed vector that is equal or is less (i.e., this ratio is the required p-value; for mathematical details see Li, 2011).Here, a p-value below .05indicated significant similarity in the direction of divergence trajectories.Finally, also from TestOfAngle, the angle for the divergence vector between pairwise groups was determined by calculating an arc cosine.

Results
Linear measurements in wild-caught and F1 fish Wild-caught sticklebacks from warm habitats had longer second (but not first) dorsal spines in most population pairs, as indicated by an interaction between thermal habitat and population pair (p < .001;see online supplementary material, Table 2, Figure 3).Tukey's HSD identified that this interaction was driven in part by differences that occurred between warm and cold populations within the A1 and S2 pairs (p < .05).The effect of thermal habitat on pectoral fin length also varied across population pairs (p = .005;see online supplementary material, Table 2).This also included an interaction between thermal habitat and population pair (p < .001;see online supplementary material, Table 2), whereby sticklebacks from warm habitats appeared to have similar pectoral fins in some pairs (A2 and S1) but shorter pectoral fins in other pairs (A1 and S3; Figure 3).
In contrast, lab-reared F1 sticklebacks displayed a significant effect of thermal habitat in the first dorsal spine length but not in the second spine (see online supplementary material, Table 2).Lab-reared F1 sticklebacks from warm habitats had longer first dorsal spines relative to those from cold habitats in most population pairs (Figure 3).Our results, therefore, demonstrated general heritable differences in first dorsal spine length based on the thermal habitat, but there was no such evidence for divergence in second dorsal spine length or pectoral fin length between thermal habitats, population pair, or their interaction (see online supplementary material, Table 2).

Body shape variation across all wild-caught and F1 fish
Our nested ANOVA showed no significant difference between RRPP models that included all factors and one that was simpler by excluding allopatry/sympatry status (nested within population pair).Therefore, our best model included thermal habitat, population pair, and centroid size as factors.Thermal habitat affected body shape indicating a degree of shared divergence across all population pairs (both wild and labreared versions), but it was also apparent that interactions between thermal habitat and population pair explained substantially more morphological variation than population pair alone (Table 2).

Body shape variation within populations
Models for each of the population pairs indicated thermal habitat effects (p < .05)for all comparisons (Table 3).Notably, this included thermal habitat effects in the F1 lab-reared fish indicating that morphological differences were heritable.Effect sizes for thermal habitat, as indicated by r-squares, varied across wild populations with a range between 3% and 18% for wild fish, and between 3.0% and 15.0% in F1 lab-reared fish (Table 3).Size influenced shape variation more strongly than thermal habitat in lab-reared fish (in 3 of 4 cases) where it had significant effects.In wild fish, size was less of an influence as thermal habitat consistently had the strongest effect, and only half of the wild pairs (3 of 6) showed a significant effect of size.Divergence in allometry was also suggested by thermal habitat interactions with size in half of the wild and lab-reared comparisons (Table 3).Deformation grids indicated that wild-caught sticklebacks from warm habitats tended to have a steeper craniofacial profile, a more subterminal jaw, a larger preorbital region, and a deeper mid-body with a shorter and more posteriorly tapered caudal peduncle relative to cold fish (Figure 4).Lab-reared F1s derived from warm habitat fish tended to uphold these differences in morphology (Figure 4).

Orientation and magnitude of thermal divergence trajectories
Based on our trajectory analysis, we could not reject the null hypotheses that populations do not differ in magnitudes or orientations of divergence.However, in no case did the angle between divergence trajectories fall below 15°, thus our strict criterion for parallel evolution was not met (Table 4).In some cases (9 of 25) support for divergence vector orientation similarity was greater than expected by chance, on the basis of the closed-form approach (Table 4).Among our comparisons, the strongest support for similarity in divergence occurred between the lab-reared versions of S1, S2, and A2 population pairs where all pairwise comparisons showed support from the closed-form approach.Also, for these pairs, the null was not rejected from our trajectory analysis, while the two smallest angles between trajectories overall (55.0°, 60.7°) were found among these comparisons (Table 4) Across wild-caught fish, our trajectory analysis tests for differences in pairwise path distances of thermal habitat divergence trajectories revealed five cases where magnitudes of divergence differed (Table 4).This included three comparisons Figure 3. Box and whisker plots depicting variation in linear traits with significant effects from thermal habitat and population pair in wild and F1 labreared sticklebacks.Each of the ambient (blue) and geothermal (red) populations are represented with allopatric (designated by an "a") and sympatric (designated with an "s") groups named following the details provided in Table 1.Linear traits were regressed against centroid size to generate residuals, with each box representing the median, upper, and lower quartiles.Panel (A) depicts variation in pectoral length in wild fish, (B) variation in spine 2 length, and (C) variation in spine 1 length within lab-reared fish.Evolution (2023), Vol. 77, No. 1 between wild population pairs, but where pairs differed substantially in age.The magnitude of divergence was larger in the younger rather than the older pair [A1 (0.044 distance) vs. A2 (0.015 distance), A1 vs. S3 (0.023 distance), and A2 vs. S2 (0.036 distance)] in all cases.The remaining two cases of divergence in magnitude involved wild versus lab comparisons where in both cases divergence in wild fish exceeded that found in lab-reared fish (A1 wild = 0.044 distance, A1 lab = 0.014 distance, S2 wild = 0.036 distance, and S2 lab= 0.020 distance).

Discussion
We tested for morphological divergence between sticklebacks from different thermal habitats and examined whether such divergence followed shared patterns across six population pairs.We found consistent evidence for divergence in relation to thermal habitat, but shared differences in body shape (while significant) appeared to have a weaker influence relative to population effects and their interaction with thermal habitat.However, based on our pairwise trajectory analysis, some population pairs from our common garden experiment showed evidence for shared effects from thermal habitat divergence.In these instances, we observed that angles between population pairs were smaller in the lab-reared fish relative to their wild counterparts.This can suggest that shared patterns of plasticity could be involved, in that fish from the two respective thermal habitats respond similarly to lab conditions.This could suggest that plasticity has a common ancestral basis across these populations, but it could also be that divergence in plasticity (in response to lab and wild conditions) between cold and warm populations has occurred in a similar way across population pairs.In this study, plasticity is revealed by differences between lab and wild fish, but this means attributing such differences to a single factor is not possible given the numerous ways these environments can differ.Further experiments that focus on plasticity across different rearing temperatures could be especially useful for addressing these different possibilities, especially for testing whether genotype-by-environment interactions have evolved between thermal habitats.
We also found that connectivity between thermal habitats had no significant role in thermal habitat divergence within our model that included all study fish.Nonetheless, allopatric pairs did display smaller angles between divergence trajectories relative to sympatric population pairs and some congruence between tests of similarity in divergence (Table 4).Thus, while not conclusive, our findings are indicative that allopatry may result in a greater degree of predictability in response to warmer habitats.

Effects of thermal environment on morphological evolution
It is well established that water temperature can influence body shape development in fishes through within-generation plastic responses (Marcil et al., 2006;Ramler et al., 2014;Sfakianakis et al., 2011).Generally, higher temperatures lead to increased body depth (Marcil et al., 2006;Rowiński et al., 2015;Sfakianakis et al., 2011).Consistent with this, we found that in most population pairs, wild-caught sticklebacks from warm habitats were deeper-bodied than those from cold habitats.They also tended to have steeper craniofacial profiles, a more subterminal mouth, and a longer second dorsal spine (Figures 3 and 4).
Although our study presents evidence for heritable morphological differences that appear stereotypical for these traits between sticklebacks from different thermal habitats, it is unclear whether these differences are adaptive.Since the effects of temperature on body shape can be either direct developmental responses or indirect (i.e., mediated by changes in other ecological conditions), the observed morphological divergence between thermal habitats could be due to selection from such indirect effects.For example, changes in jaw orientation and body depth may be driven by differences in food availability (Rowiński et al., 2015) or diet composition (Hjelm et al., 2001).More specifically, greater body depth and steeper craniofacial profiles can be related to swimming and jaw function, as it signifies hypertrophied epaxial musculature.This suggests an increased suction ability and, along with a subterminal jaw, would indicate a more benthic foraging ecology in the warm-habitat sticklebacks (McGee et al., 2013).Finally, it is also possible that changes are occurring through nongenetic sources of inheritance.Indeed, transgenerational effects have been documented for traits related to life history, behavior, and physiology in this species (Afseth et al., 2022;Shama et al., 2014), making it highly plausible that other traits such as morphology could be similarly affected.
Another potential explanation for our findings is that sticklebacks in warm habitats have evolved deeper bodies in response to higher predation risk.Greater body depth is thought to improve predator escape performance increased maneuverability or predator gape limitations (Domenici et al., 2008;Reimchen, 1991;Walker, 1997).
Similarly, dorsal spines are an antipredator defense and are generally longer in populations that experience elevated predation pressure (e.g., Blouw & Hagen, 1984).Here, we found evidence for longer second dorsal spines in wild fish from warm habitats but longer first dorsal spines in lab-reared F1 fish from warm habitats.This suggests that evolved differences in the developmental systems between thermal habitats compensate to keep the first spine the same under variable natural conditions but allow the second dorsal spine to respond via plasticity.In our study system, we expect bird predation to be higher in warm habitats due to the lack of ice cover during the winter and the fact that birds tend to be attracted to warmer areas (Rowiński et al., 2015).On the other hand, sticklebacks in warm habitats likely experience a lower risk of predation from freshwater piscivorous fish, which may be unable to cope with high temperatures (Eliason et al., 2011).Further research will be needed to investigate the

Magnitude of thermal divergence in body shape
The magnitude of divergence between populations in warm and cold habitats may be influenced by various factors.For example, some researchers have argued that gene flow will constrain divergence between ecotypes when there is a potential for dispersal (Hendry & Taylor, 2004;Lenormand, 2002;Slatkin, 1985).Under this scenario, we would expect sympatric population pairs to be less divergent than allopatric population pairs (Hendry & Taylor, 2004;Pinho & Hey, 2010).However, sympatric pairs could instead be more divergent because of character displacement, whereby differences between ecotypes are more pronounced in areas where they co-occur and minimized in areas where their distributions do not overlap.This pattern results from trait evolution driven by competition among ecotypes, or closely related species, for a limited resource (Brown & Wilson, 1956;Losos, 2011).
Here, we found three cases where the magnitude of divergence differed between population pairs and two of these involved larger magnitudes of divergence being present in allopatric pairs relative to sympatric pairs.Nonetheless, neither the presence or absence of geographical barriers nor population age seemed to consistently influence the magnitude of thermal habitat divergence in the body shape of wild-caught sticklebacks (Tables 3 and 4).This suggested that the degree of response to thermal habitat variation may be population and site-specific rather than generalizable against these factors.If so, it may mean that degree of morphological response to climate change will be difficult to predict.
Nonetheless, it could be expected that for populations that have been diverging for longer, there would be more scope for natural selection and genetic drift to introduce adaptive or stochastic phenotypic differences (Ord & Summers, 2015).We, therefore, expected relatively young population pairs (<100 years old) to be less divergent than old population pairs (>1,000 years old).However, this was not the case, with population pairs that have existed for just a few decades being equally or even significantly more divergent than population pairs that have existed for thousands of years (Tables 3  and 4).This suggests that the main aspects of thermal habitat divergence could evolve rapidly and reach an "evolutionary plateau," whereby finer scale change occurs at a slower rate.Indeed, across a range of taxa, it has been shown that contemporary evolution initiates with an extremely rapid rate of change that scales negatively with time (Stockwell et al., 2003).This is also the case with other examples of parallel morphological evolution in fish, such as in adaptive radiations of African cichlids (Cooper et al., 2010).Our findings suggest that evolutionary divergence, while not strictly parallel, is somewhat consistent different timescales.
Mechanistically, this could be in line with the predictions of Orr (2005) whereby a small number of beneficial mutations within a common ancestor enhances the probability of parallel evolution under a common selective environment.Much remains to be explored about the genetic relationships and mechanisms involved with the evolutionary divergence we have identified, as well as its adaptive value, but our findings at least indicate that fishes could rapidly evolve in response to climate change.

Direction of thermal divergence in body shape
Habitat connectivity may influence not only the magnitude of divergence but also the degree of parallelism between replicate populations (Bolnick et al., 2018).Gene flow between different habitat types is thought to constrain local adaptation within each habitat, so if there is variation in the extent of gene flow among replicate populations, migration-selection balance will act differently contributing to nonparallel evolution (Hendry & Taylor, 2004;Moore et al., 2007;Stuart et al., 2017).As a result, allopatric population pairs (i.e., with physical barriers to gene flow) and sympatric population pairs (i.e., with more potential for gene flow) may differ in their degree of parallelism.Indeed, we found some evidence that geographical isolation may facilitate similar trajectory orientations in warm-cold population pairs of wild-caught sticklebacks: two of three allopatric population pairs had evidence of similar trajectories across methods, whereas no sympatric population pair comparisons showed congruence across methods (Table 4).Notably, congruent evidence across methods did not persist in F1 lab-reared versions of these allopatric pairs suggesting that development for shared patterns of divergence partially relies on natural environmental cues.Previous theoretical and empirical work suggests that the degree of parallelism between replicate populations may also depend on the duration of evolutionary divergence (Lucek et al., 2014;Ord & Summers, 2015).For example, population pairs that have been diverging for longer have more scope for natural selection and genetic drift to alter their evolutionary trajectories, resulting in a lower degree of parallelism in older populations (Bolnick et al., 2018).Yet, if evolution is limited by mutation rate, older populations will have had more time to accumulate similar adaptive mutations that produce a similar phenotypic solution in response to a particular environment (Orr, 2005;Whitlock & Gomulkiewicz, 2005).Under this scenario, older populations would have a higher degree of parallelism.Our findings did not support either of these possible outcomes, since the extent of parallel divergence did not differ consistently between young and old population pairs.Given the overall young age of all the populations, we investigated the accumulation of lineage-specific mutations that are unlikely to play a role.However, it could simply be that divergence is driven primarily by specific local conditions as evidenced by the larger effects of population pair and their interaction with thermal habitat relative to thermal habitat alone.
It should be noted that even though we have discussed the effects of habitat connectivity and population age, several other factors could influence the magnitude and direction of divergence in wild populations.These include ancestry and evolutionary history (Langerhans & DeWitt, 2004), initial and ongoing effective population sizes that could cause differential bottlenecks across populations (Szendro et al., 2013), variation in sexual selection (Bonduriansky, 2011;Maan & Seehausen, 2011), and many-to-one mapping, which refers to multiple phenotypic solutions to the same functional problem (Gould & Lewontin, 1979;Wainwright et al., 2005).Given our evidence for broad qualitative similarities in divergence (Figure 4), as well as our quantitative results, further work is needed to examine the contribution of these additional factors.

Predictability of evolution and adaptation to climate change
As discussed above, we found evidence for thermal habitat divergence between all population pairs.While strict parallelism was not supported, there is evidence for some degree of consistency across population pairs, which could be attributed to natural selection, developmental bias, or their interaction (Brakefield, 2006;Losos, 2011;Uller et al., 2018).We cannot separate the effects of these processes in the present study, but regardless of the underlying causes, our results suggest that morphological responses to increasing temperatures could be partially predictable for fish populations.Under these conditions, we may expect fish to evolve subterminal jaws and a deeper mid-body after being exposed to elevated temperatures over multiple generations, a noteworthy trend in the evolution of these populations given similar observations from previous studies on morphological plasticity in fish (e.g., Marcil et al., 2006;Sfakianakis et al., 2011).On the other hand, the migration of individuals between different thermal habitats or microhabitats may exaggerate nonparallel evolution (Bolnick et al., 2018;Oke et al., 2017) and reduce our ability to predict evolutionary responses to climate change.
At the same time, our findings suggest that a large degree of divergence in response to thermal habitats is driven by nonparallel, population-specific divergence.It is possible these aspects of morphological change are important for local adaptation, as the natural ponds and lakes containing our focal populations were not experimental replicates and varied in size, depth, and likely many other biological aspects.Also, these unique aspects could be driven by local population history, but it appears that all of the populations we used were derived from a common ancestor (Costa et al., unpublished SNP data).Future work should aim to address whether both shared and unique aspects of divergence have adaptive impacts.While both could contribute to adaptation, there may be a sequence of events whereby adaptive variation driven by developmental bias, which initially enhances parallelism (Uller et al., 2018), is accessed more readily in the early stages of divergence.Evidence is emerging to suggest that plasticity could provide this type of initial evolutionary response (Noble et al., 2019;Parsons et al., 2011) with specific fine tuning (i.e., local adaptation) occurring at later stages.This appears to be the case in other examples of adaptive radiation (Cooper et al., 2010).However, given that our findings do not suggest population age is associated with the degree of similarity in the orientation of divergence, it could be that local adaptation is able to happen in concert with parallelism.To address this, future work could examine the range of variation available in ancestral cold and marine populations, especially with respect to how they respond to changes in temperature.

Figure 1 .
Figure 1.Map of Iceland showing the sampling locations of warm-and cold-habitat sticklebacks we collected for this study.All sticklebacks were collected from freshwater populations, and each of the six population pairs (A1, A2, A3, S1, S2, and S3) is indicated by a different color.

Figure 4 .
Figure 4. Deformation grids depicting shape differences (3× magnification) between the average geothermal (red) and ambient (blue) fish within population pairs under natural wild, and lab conditions.Allopatric population pairs are designated as A1, A2, and A3, while sympatric pairs are S1, S2, and S3.Outlines connect the outermost landmarks, and landmarks surrounding the eye for each group.

Table 1 .
Sampling locations and sample sizes of warm-and cold-habitat sticklebacks collected in May and June of 2016.

Table 2 .
Results of an RRPP model testing the effects of thermal habitat, population pair, and their interaction on stickleback shape from wild-caught and F1 lab-reared sources.Df denotes degrees of freedom, SS denotes sums-of-square, and Rsq denotes the r-square for each factor.Statistically significant p-values are indicated in bold.

Table 3 .
Results of separate RRPP models testing the effect of thermal habitat and geometric centroid size on each population pair in wild-caught and F1 lab-reared sticklebacks.
functional significance and the developmental mechanisms (i.e., what allows for or causes canalization) of the morphological differences we have documented.

Table 4 .
Trajectory analysis results from an RRPP model on the body shape of wild-caught and F1 sticklebacks from cold and warm thermal habitats.represents an additional test of nonorthogonal trajectory based on closed-form formulas at a critical angle below 75.6° and indicates significant similarity in divergence trajectory orientations when the alpha value is below 0.05.Magnitude differences in path distances for each pair are also provided and p-values for pairwise tests of path distances are provided in the far-right hand column.Statistically significant values are indicated in bold.
Here R denotes slope vector correlations for each pairwise comparison.p-values are provided for the RRPP-based tests of trajectory similarity.Additionally, for trajectory orientation P(HC)