{'fraction_non_alphanumeric': 0.02639253893981155, 'fraction_numerical': 0.014550349336580988, 'mean_word_length': 5.250500801282051, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 1, 'zh_pornsignals.json': 0}}
Dimensionality and modularity of adaptive variation: Divergence in threespine stickleback from diverse environments

In nature, populations are subjected to a wide variety of environmental conditions that affect fitness and induce adaptive or plastic responses in traits, resulting in phenotypic divergence between populations. The dimensionality of that divergence, however, remains contentious. At the extremes, some contend that populations diverge along a single axis of trait covariance with greatest availability of heritable variation, even if this does not lead a population directly to its fitness optimum. Those at the other extreme argue that selection can push populations towards their fitness optima along multiple phenotype axes simultaneously, resulting in divergence in numerous dimensions. Here, we address this debate using populations of threespine stickleback (Gasterosteus aculeatus) in the Cook Inlet region of southern Alaska from lakes with contrasting ecological conditions. We calculated effective dimensionality of divergence in several trait suites (defensive, swimming, and trophic) thought to be under correlated selection pressures, as well as across all traits. We also tested for integration among the trait suites and between each trait suite and the environment. We found that populations in the Cook Inlet radiation exhibit dimensionality of phenotype high enough to preclude a single axis of divergence.


Introduction
Evolutionary biologists have long grappled with uncertainty regarding the dimensionality of evolutionary processes and patterns, which has important implications for understanding rates and patterns of evolution, as well as their constraints. Although the dimensionality of evolution can be considered from selection oriented or genetic variance oriented perspectives, it can be broadly simplified to two extremes. At one extreme, divergence is constrained by the major axis of the genetic covariance matrix (the "G matrix"), to a "genetic line of least resistance" (Gmax), irrespective of the strength, complexity, or direction of selection. At the other extreme, highdimensional divergence is expected to be facilitated both by a complex selection landscape influencing mostly uncorrelated traits with genetic variance distributed across them (Mezey and Houle 2005;Nosil et al. 2009). (Schluter 1996;Hine et al. 2014).
The dimensionality of evolution has important implications for rates of evolutionary change, the directions in trait space in which adaptation is possible, and the number of niche spaces that can be filled. For instance, if the dimensionality of traits is low and correlated traits experience opposing selection pressures, the rate of evolution will be lower than if selection was focused in a single direction similar to the main axis of genetic variation (Agrawal and Stinchcombe 2009). By contrast, high dimensionality of traits could mean fewer constraints on the number of possible combinations of trait values, and thus a greater number of ecological niches that can be rapidly filled by adaptation (Kirkpatrick 2009). Further, multidimensionality of diversification also promotes many-to-one mapping of phenotype to function; that is, the existence of multiple possible phenotypic "solutions" to selection that results in similar fitness Wainwright et al. 2005), and accommodates the ability of populations to adapt to different peaks on a rugged fitness landscape (Langerhans 2018). In this sense, highdimensional phenotypic divergence is more aligned with the empirical evidence on the limits of parallel evolution to similar selection pressures than unidimensional divergence (Oke et al. 2017;Langerhans 2018).
One way of conceptualizing and quantifying evolutionary dimensionality is through the related concepts of modularity and integration. Modules are defined as groups of traits that are strongly correlated with each other, but less strongly correlated with other groups of traits (Wagner and Altenberg 1996;Mitteroecker and Bookstein 2007;Wagner et al. 2007). Integration describes correlations between traits or between trait groups that do not permit independent variation (Cheverud 1982;Wagner and Altenberg 1996;Mitteroecker and Bookstein 2007;Adams and Collyer 2016). Integration between modules can constrain evolution by limiting the number of possible evolutionary trajectories without reorganization of the G-matrix (Arnold 1992). A high degree of modularity allows for more flexibility in evolutionary trajectories because of the relative independence of each trait group (Wagner and Altenberg 1996), although it has been predicted that dimensionally complex evolution may proceed at a slower rate than evolution along a single axis (Welch and Waxman 2003). However, measuring dimensionality of evolutionary change is made more difficult by patterns of modularity, which results in its uneven distribution across traits.
Most previous studies have tended to focus on associations within one or two trait modules of functionally similar traits likely experiencing similar selection pressures, eliminating the possibility of exploring how integration and modularity expand or contract the potential adaptive range of the phenotypes. For example, Schluter (1996) focused on traits related to trophic ecology in sticklebacks, Darwin's finches, and sparrows. Hine et al. (2014) considered only cuticular hydrocarbons in Drosophila. Hansen & Houle (2008) used only different measurements of Drosophila wing shape. Kirkpatrick (2009) usedin separate analysesfat and muscle traits for beef cattle, body shape measurements in the fish M. eachamensis, wing shape measurements in Drosophila, lactation curves for dairy cows, and growth trajectories for mice. We suggest that these previous assessments of associations among functionally similar traits are perhaps more likely to represent one or only a few trait modules, which would result in an underestimate of evolutionary dimensionality and an overestimate of the constraints on evolution. Hence, our goal in the present study is to assess the extent of the dimensionality of evolutionary divergence within and between hypothesized modules (hereafter "trait suites") with different selection pressures, and the degree to which the measurement of dimensionality is influenced by trait modularity and integration. We did this by considering multiple traits in each of three trait suites of different functional categories in evolutionarily divergent populations, making it unlikely that apparent restriction of phenotypic divergence primarily to a single dimension is a consequence of the traits being under similar selection pressures.

Study system and questions
Threespine stickleback (Gasterosteus aculeatus) are an excellent species for studying evolutionary dimensionality through the lens of modularity and integration because they are highly variable within and among populations, and because multiple freshwater populations are independently derived from ancestral marine populations (Bell and Foster 1994). Many of their variable traits are known to be correlated with environmental variables (Bell et al. 1993;Reimchen 1994;MacColl et al. 2013;McGee et al. 2013) that can be used to classify traits into functional suites based on the source(s) of selection they experience. Furthermore, stickleback divergence between habitat types (e.g., lake-stream (Lavin and McPhail 1993), benthic-limnetic (McPhail 1993), and marine-freshwater (Gelmond et al. 2009)) has been described as occurring along axes corresponding to habitat, but whether a single axis can accurately describe any of these divergence patterns is uncertain.
Here, we ask how much variation in phenotype can be explained by a single axis of divergence, and whether phenotypic variation among variable habitats is better explained by a multifarious or unidimensional model of trait divergence. If trait divergence is close to unidimensional, then it may be constrained either by similarly low-dimensional selection, or by strong genetic or developmental integration between traits. Our study focuses on lacustrine environments, so the benthic-limnetic axis of divergence is likely to represent a major axis of phenotypic differentiation, and will occasionally be referred to give an impression of the differences being described. This axis can be defined by a stouter and deeper body shape, a shortened head with a downturned jaw, and fewer gill rakers at the benthic end, and more elongate, zooplanktivory-adapted fish at the limnetic end (Bentzen et al. 1984;McPhail 1993;Willacker et al. 2010). If divergence is multidimensional, it will occur along multiple axes that are selected upon by different environmental variables.
We collected stickleback and conducted environmental surveys at 14 lakes in the Cook Inlet region of southern Alaska. We divided trait measurements of the collected sticklebacka prioriinto four functional suites (defensive, trophic, swimming, body shape) that could be under differing selection pressures or respond plastically to different environmental conditions. We then used these morphological data to answer four main questions ( Figure 1).

1)
To what extent do lacustrine stickleback diversify along a single axis of divergence among populations, and what is the effective dimensionality of that divergence? Although there is no formal null expectation of effective dimensionality for this question, we predicted that populations of stickleback confined to lakes diversified along a single axis of correlated variation between phenotypes, and that divergence is largely a unidimensional process. We applied this question to each of the pre-defined functional trait suites, as well as a single dataset including all traits. If there is a single primary axis of divergence, the primary axes of the trait suites will be correlated. If trait suites are uncorrelated, then dimensionality of divergence is not constrained to the major axis of the G matrix unless the trait suites themselves are highly multidimensional.
2) To what extent do lakes share common within-lake correlation structures between trait modules? We predict that if trait suites are closely correlated, then lakes will share common within-lake correlation structures between trait suites, reflecting conserved covariance among traits suites.
3) To what extent do functional trait suites covary as adaptive modules? Trait suites that constitute modules would strongly covary internally, and be representative of axes of divergence that could be predicted based on the environmental conditions imposing selection on the module, irrespective of pressures on other traits or trait suites.

4)
To what extent are different environmental variables, or groups of environmental variables, associated with variation in all phenoptypes and in the trait suites? Phenotypic divergence corresponding with environmental variation would indicate that the phenotypic variation is likely adaptive or shares common plastic responses, and that it is only weakly constrained by correlations between trait suites.
From previous work, we do have some a priori expectations of trait responses to the environment. For defensive traits, low dissolved calcium concentration in the absence of predators is strongly predictive of reduction of the pelvis and spines (Giles 1983;Bell et al. 1993). Although some of the lakes we sampled have stocked or native salmonids that could act as predators, only limited information about the predator communities of most of the lakes is available (ADF&G 2015). Robustness of defensive traits is also associated with larger, clearwater lakes with higher visibility for predators (Reimchen et al. 2013), and these lakes tend to have fish with narrower caudal peduncles and body shapes with more anteriorly-shifted dorsal spines and dorsal and anal finstraits that are not directly related to defense (Spoljaric and Reimchen 2007). Willacker et al. (2010) also found that depth-related variables were significantly correlated with the benthic-limnetic axis for skull morphology. One populationthat of G Lakewas known to have an introduced population of Muskellunge (Esox masquinongy; Lipka et al. 2020;Massengill et al. 2020), a congener of Northern Pike (E. lucius), which Willacker et al. (2010) found was associated with more limnetic cranial morphologies in stickleback than those in populations from otherwise comparable habitats. We expect that lakes with environmental variables that depart from the mean lake in more dimensions will have stickleback that are divergent in more directions. We expect that trophic traits could depend on Among-lake divergence is uncorrelated among traits suites Single major axis of among-lake trait divergence. Trait divergence is lowdimensional.
Multiple axes of among-lake trait divergence. Trait divergence is highdimensional.
Among-lake divergence is correlated between trait suites Trait suites are highly integrated, with common integration structures between populations Trait suites are highly modular (evolve independently) Trait suites are integrated within lakes, but correlation structures are not shared among lakes.

Variation between environments is associated with variation in traits between populations
The relationship between environment and phenotype differs by trait suite and along different ecological dimensions Environmental variation is not associated with variation of traits Q1b -Are the patterns of divergence correlated between trait suites? Q4 -Is environmental variation between populations associated with phenotypic divergence?
Q2&3 -To what extent are trait suites modular or integrated with each other, and do integration patterns differ between populations?
Q1a -Is dimensionality of divergence high or low among all traits? This question is also applied to trait suites. Figure 1. Conceptual diagrams illustrating the potential results of the questions addressed in this paper. Grey ellipses represent total variation in phenotype between populations. Filled, colored ellipses represent variation within trait suites between populations, open dashed ellipses represent populations, and solid circles represent population means. Lines bisecting population ellipses indicate the major axis of variation within that population for the trait suite in which it is shown. Panel axes, except for in D1-D3, are arbitrary dimensions in a shared trait space, but trait suite ellipses are separated in some panels for illustrative purposes. For the tests that correspond with the panels, please see Table A1.
Phen Phen nutrient availability, mediated by the prey community. However, because of the dynamic feedbacks between prey availability and trophic morphology, the direction of causality between morphology and some environmental variables that influence predator-prey eco-evo relationships is not clear (Brooks and Dodson 1965;Post 2008, 2009;Schmid et al. 2019).

Study system
In This design increased the likelihood that populations had sufficient phenotypic diversity to detect multidimensionality of divergence if it is present in wild populations more generally, and although some of these populations are close to one another geographically, they do not share surface connections that are plausible routes of dispersion. Generally, in the Cook Inlet region of Alaska (Willacker et al. 2010)as elsewhere in the species range (Schluter and McPhail 1992), "benthic" stickleback are adapted to feeding on benthic macroinvertebrates and have shorter heads and deeper bodies, whereas "limnetic" stickleback are adapted to feeding on zooplankton, and have more elongate body and head shapes. for analysis). Using the digital photos, the following measurements were made in ImageJ (Rasband 1997(Rasband -2018: caudal peduncle width (CP), upper jaw length (JL), snout length (SN), head length (HL), and eye diameter (ED) (Figure 3).
Linear measurements were adjusted for allometry using the formula Mstd = MO ( ̅̅̅ / SLO) b , where MO and SLO are the observed trait length and standard length, respectively, ̅̅̅ is the  (Reist 1986;Hendry and Taylor 2004). GR, LP, and LPM were not adjusted for allometry because all but 2 of our 792 fish had standard lengths < 30 mm, the approximate size by which all gill rakers, lateral plates, and pelvic bones have formed (Bell 1981(Bell , 1987(Bell , 2001Glazer et al. 2014). Spearman's rank correlations, with p-values adjusted for multiple comparisons using the false detection rate (FDR) method (Benjamini and Hochberg 1995), confirmed that GR and LP were not correlated with SL in any lake.
Following Stuart et al. (2017), univariate measurements were then assigned to one of three trait suites (Defense, Trophic, Swimming) based on the type of selection they are expected to experience most directly (Table S1). Although some traits could be associated with more than one trait suite based on selection in the environment, each had to be assigned to a single trait suite for methodological reasons. Additionally, the tests of correlation between trait suites described below implicitly acknowledge the potential for functional overlap between trait suites.
Missing values of traits, except GR, which had by far the lowest N of all traits recorded (NGR = 276, or 34% of specimens), were estimated using the regression formula of the most correlated trait with r > 0.1 within the same trait suite (numbers of individuals for which traits were thus estimated are provided in Table S2). Units of traits were then standardized for multivariate analyses by transforming units into standard deviations of the total sample, and centered about zero (Z-transformation). All data analyses were performed using the R statistical environment, version 3.6 (R Core Team 2019).

Geometric morphometrics
Body shape was quantified using landmark-based geometric morphometric (GM) methods, with thirteen homologous anatomical landmarks placed from mouth to tail, to which eight semi-landmarks were added to capture the curvature of the head from the tip of the snout to the first dorsal spine (Figure 3). Landmarks and semi-landmarks were placed using tpsDig2 2.31 (Rohlf 2018). A Procrustes superimposition was performed on all specimens using the gpagen function in the R package geomorph , with curves defined using minimum bending energy of the semi-landmarks. Photographs of Long Lake specimens taken with the two different cameras were found to be significantly different with a Procrustes ANOVA, so landmark data of specimens photographed with the Nikon D800 was excluded from geometric morphometric analysis. Specimens for which the location of the landmark at the base of the first dorsal spine was uncertain were excluded, as this landmark was used to place the series of eight semi-landmarks. Repeatability of shape was assessed to be 92% using variance ratios from Procrustes ANOVA as described by Zelditch et al (2012).
Because the GM data was available for two fewer lakes than the data describing the other trait suites, the subsequent analyses were conducted both including and excluding the GM data, with results provided in the Supplementary Tables when they are not provided in the Results section. Although shape is described as a trait suite throughout this paper, it is expected to be highly integrated with all three of the other trait suites. This is because it is associated with Defense from gape-limited predators (Reimchen 1991), with Defense and Swimming traits via its association with the fast-start response for predator avoidance (Taylor and McPhail 1986;Bergstrom 2002), and with Swimming and Trophic traits via the relationship between body depth and maneuverability in spatially complex foraging environments (Walker 1997;Feilich 2016).
Additionally, linear measurements included in other trait suites, especially BD, ED, CP, and HL describe aspects of shape, and thus these morphological traits are inextricably linked with the shape data.

Limnological Conditions and Dietary Indicators
Lake

Question 1 -Dimensionality of divergence, and correlations of trait suites
We examined the extent to which trait divergence among lake populations was multidimensional by first applying a linear discriminant analysis (LDA) including traits from all trait suites Ripley et al. 2020 However, because these categories ease discussion, we regard "low" effective dimensionality as < 2, "moderate" from 2 -3, and "high" > 3. These values were chosen because the dimensionality of a circle projected onto a plane is 2, and the dimensionality of a sphere is 3, making their interpretation somewhat more intuitive than other thresholds. For convenience, we will refer mostly to the Dn1 effective dimensionality index unless otherwise noted. Because we are testing the extent to which this radiation aligns with to hypotheses along a continuum of dimensionality, there is no null hypothesis, and expectations of dimensionality are likely to vary by system and environmental complexity.
Because gill rakers were counted in only approximately one third of specimens, the above discriminant analysis for all traits was conducted with and without gill rakers. Additionally, this discriminant analyses excluded geometric morphometric landmark coordinates, both because landmarks excluded two lakes and because some components of shape are redundant with other variables.
As with the dataset of all traits, we conducted LDAs within Defense, Swimming, and The effective dimensionality is limited by the number of traits included in each trait suite (e.g., the defensive trait suite, which includes 5 traits, and cannot have an effective dimensionality greater than 5), and so we implemented a rarefaction-like permutation procedure to visualize changes in effective dimensionality of trait suites as traits are accumulated at random, and compare suites of different numbers of traits. For details on this procedure and results, see Appendices A and B, respectively.

Question 2 -Comparisons of within-lake correlation structures
We used two-block partial least squares correlation (2B-PLS), to test whether the trait suites were correlated within populations. This technique takes two data matrices from the same individuals, for example our data for swimming and defensive trait suites, and identifies the linear vectors through each dataset that maximize covariation explained between the datasets (Rohlf and Corti 2000). We conducted 2B-PLS tests for all pairwise comparisons of trait suites.
Each test yielded both a correlation coefficient (rPLS), and a Z-score ). The Zscore serves as an effect size of integration between the two blocks (here, trait suites) for each population, and which is comparable between 2B-PLS tests on populations with different sample sizes (Adams and Collyer 2016). We used these standardized effect sizes to compare integration strength of trait suites between populations with pairwise two-sample Z-tests (Adams and Collyer 2016;, and characterized the orientation of those differences by calculating the angle between vectors describing the most variance through data blocks between lakes. 2B-PLS tests and between-population Z-tests were performed using the two.b.pls and compare.pls functions, respectively, in the R package geomorph ). Reported p-values were adjusted for multiple comparisons using the FDR method (Benjamini and Hochberg 1995).
By way of a simple example, in a radiation of three populations (A, B, and C), we could conduct 2B-PLS tests between defensive traits and swimming traits in each population. The rPLS statistic would give us the level of correlation between the defensive traits and swimming traits.
Because rPLS is sensitive to sample size, however, we compare the strength of correlations between trait suites using Z-scores which describe the strength of the effect relative to the null expectation of no association. If, after conducting pairwise Z-tests between populations, the comparisons between populations A and B yield a low effect size, but the A-C and B-C comparisons yield high effect sizes, then the strength of integration between defensive and swimming traits is similar between A and B but C differs from both. In this example radiation, then, the strength of the relationship between defense-associated and swimming-associated traits in population C diverged from that of the other two populations. Calculating vector angles between data blocks of the same trait suite for pairwise comparisons between lakes allowed us to compare orientation of integration structures, which are not taken into account in Z-tests for differences in strength of effects between lakes.
Because PLS methods are sensitive to both vector magnitude and orientation, our use of correlation matrices does produce different results in our 2B-PLS tests than use of covariance matrices would, (and indeed, we did use unstandardized GM landmarks in these analyses, to the same effect). However, use of a correlation matrix for instances in which data do not all exist on the same scale, as is the case here, is recommended by Rohlf and Corti (2000), and even required for uses of 2B-PLS in which a block of morphological data is associated with a block of environmental data or other non-morphological data related to niche occupation, as environmental data are rarely on the same scale (Corti et al. 1996;Fadda and Corti 1998;Felice et al. 2019).

Question 3 -Evolutionary modularity of trait suites
We tested the extent of modularity of trait suites, both within each population and across all populations using the modularity.test function of geomorph, which tests whether the ratio of covariation between to within designated trait suites is less than if traits had been assigned to suites randomly (Adams 2016;. When this test is applied to multiple trait suites simultaneously, the average of pairwise covariance ratios (CR) is used as a test statistic (Adams and Collyer 2019), which provides a better picture of overall modularity, but provides less information about the independence of each trait suite. Consequently, these tests were executed in two ways, involving how partitions between trait suites were set. First, we performed the test after defining a partition separating a target trait suite from all other traits.
These tests are referred to in results as "focused" tests. We then performed a test of modularity with all trait suites partitioned from one another, which we refer to in results as a "complete" test.
Because the measurements in the swimming trait suite were important components of shape, modularity tests were performed both with and without the inclusion of shape as a module to be sure that this redundancy did not alter results. All modularity tests were performed with 5000 permutations. The p-values reported for these tests are FDR-adjusted for multiple comparisons.

Question 4 -Responses to the environment
We quantified environmental differences between lakes using principal components analysis (PCA). This was first performed first using physico-chemical variables -pH, DOC, Ca, TN, TP, Maxium Depth, chl-a, and surface areafor which all data were available for all but two lakes (Echo and Ruth). This was followed by a PCA including both physico-chemical variables and stickleback foraging-related characteristics of the lakes, for which only a subset of lakes were surveyed. Both environmental PCAs used scaled environmental variables, and thus were calculated using correlation matrices rather than covariance matrices. Stickleback foraging- We subsequently executed 2B-PLS tests using normalized population-mean morphological variables as one block, and environmental variables as the other. As with PCAs, we performed 2B-PLS tests twice, using both physico-chemical and full environmental datasets, the first including more lakes (n = 12), but the second included data related to foraging (n = 9).
2B-PLS is a balanced analysis which does not assume a direction of causality. It therefore may be a more accurate characterization of our system, in which eco-evo feedbacks potentially play a role, than another multivariate regression technique like redundancy analysis (RDA), which includes multivariate regression as its first step (Legendre et al. 2011

Axes of divergence for all traits and within trait suites
In the linear discriminant analysis including all defensive, trophic, and swimming traits, dimensionality was calculated to be high ( 1 = 3.26 and 2 = 2.01) when gill rakers were included (N = 241; Tables S3, 2). In this LDA, the primary axis explained 69.0% of the variance between populations, and the first three axis explained 87.2% of the variance. The traits onto which the first axis were most heavily loaded were standard length and pelvic spine length (  An LDA of defensive traits showed that variation between populations was driven mostly by a single axis that accounted for 92.1% of this variation. Pelvic spine length (PSL) was the variable most heavily loaded onto this axis, and had a loading nearly 8 times higher than DS2L, the trait with the next highest loading. this result seems to be driven by the G Lake population, which was entirely lacking pelvic spines, and Echo Lake, which had only a single fish with pelvic spines. Dorsal spine lengths (DS1L and DS2L) drove most of the remaining variance in defensive traits between populations, and dimensionality of defensive trait divergence was low (Tables 2 and S5 Trophic trait divergence was much more multidimensional than swimming or defensive trait divergence, with effective dimensionality of 1 = 4.57 and 2 = 3.47 (Table 2). This was also much higher than De calculated from eigenvalues of individual-level or population mean trait data. In the trophic trait suite, LD1 described only 44.5% of variation between populations with PW, followed closely by HL and SN, being the most heavily loaded variable onto this axis (Table S7).
Body shape was the most multidimensional of the trait suites, with high effective dimensionality of 1 = 7.16 and 2 = 5.27, albeit calculated from data with a much higher number of traits, with each landmark having an X and Y coordinate (Table 2). For shape, the first discriminant axis only described 35.2% of variation between populations, and five LD axes were required to explain at least 80%. The first axis here describes mostly body depth, but also angle of the head and the posterior extent of the dorsal fin. As with the swimming and to some extent the trophic suite, LD1 appears to reflect the benthic-limnetic continuum ( Figures S1 and S2). Notes: The final group of values represents all traits excluding shape, because a) two populations were excluded from geometric morphometric analysis and b) our swimming traits are largely redundant with our shape data. and are shown in bold typeface to indicate that they are the De measures that are recommended for consideration, as ∞ and systematically underestimate and overestimate dimensionality, respectively.
Primary axes of all trait suite LDAs were significantly correlated with all other LD1s ( Figure S3; Table S8

Comparisons of within-lake correlation structures
Two-block PLS correlations between trait suites within populations show significant integration for each combination of trait suites in at least some lakes. Of these, the swimming x trophic suite tests show the most persistent trends (rPLS > 0.855 in all cases, Zrange = 2.36 -7.01, p ≤ 0.001 in all lakes except Engineer: p = 0.011; Table S10). The defensive trait suite and shape had by far the least consistent, and typically weakest, relationship. Their rPLS correlation coefficients ranged from 0.33 (Wik Lake, Z = -1.71, p = 0.962) to 0.79 and 0.80 (Echo Lake, Z = 3.24, p = 0.006; and Long Lake, Z = 0.99, p = 0.362, respectively). The other pairs of trait suites mostly had clear integrative relationships, although the patterns of shared variation between defensive and trophic traits were also somewhat erratic between lakes.
Pairwise comparisons between lakes for difference of PLS effect size in each pair of trait suites revealed several cases of clusters of populations that share similar integration strengths when compared with other lakes, albeit typically with relatively moderate effect sizes, and none higher than Z = 3.89 ( Figure S4). These groups do not clearly reflect populations' position along the benthic-limnetic axis, and in some cases, lakes on opposite ends of the traditional benthiclimnetic axis share similar integration strengths. Watson and Wik Lakes, for instance, have quite similar patterns of integration between defensive and trophic traits, despite Watson fish being very benthic and Wik fish being very limnetic.
More revealing were the patterns of relative orientation of PLS axes between lakes ( Figure S5). In the swimming x trophic PLS test, axes of correlation between these trait suites were strikingly parallel across all lakes, with the angles between trophic suite axes never exceeding 26º between lakes and the angles between swimming suite axes never exceeding 24º.
Integration structures between swimming and defensive traits or shape also resulted in quite parallel PLS axes for the swimming suite, although the orientations of the defensive trait especially the shape axes were much more muddled in these cases. The integration patterns between some trait suites, in contrast, sometimes clustered the lakes into more or less poorlydefined groups, within which lakes shared similar PLS axis orientations and differed in similar ways from lakes in other clusters ( Figure S5). However, the lakes composing these clusters were not necessarily the same between integration structures. Additionally and surprisingly, as was the case with integration effect sizes, the lakes in the PLS axis orientation clusters did not map clearly onto the fairly obvious, but apparently superficial, benthic-limnetic continuum. One of the clusters in the defensive block defensive x trophic suite PLS angle matrix included Wik and G Lakes, both of which are extremely limnetic, but also Tern Lake, the most benthic of any of the lakes surveyed.

Modularity of trait suites
The complete modularity test excluding shape showed that there is a statistically significant modular signal of relatively modest strength in every population except Engineer, which was only marginally insignificant, and differed only slightly from other populations (Table S11). However, the focused tests on individual trait suites revealed some patterns in modularity of trait suites. In particular, the defensive trait suite exhibited moderate to high levels of modularity in all lakes, with the effect being strongest in G, Wik, and Spirit Lakes, and weakest in Long, Finger, and Engineer (Table S12). When gill raker data is excluded from the trophic trait suite to increase sample size, the signal of modularity in defensive traits generally becomes even stronger, even exceeding Z = -10 in Wik and Spirit Lakes. Modularity in trophic traits was not consistent between lakes and exhibited a significant, moderately sized signal in four lakes (Corcoran, G, Long, and Spirit). Swimming traits were not modular in any population, and had consistently weak effect sizes. That swimming traits are not modular was not surprising given the importance of swimming performance in both feeding and avoiding predation.
The inclusion of shape in a complete modularity test lowered covariance ratios, and also slightly strengthened effect sizes of modularity (Table S13). When shape was considered in a focused test of modularity, the signal of modularity was significant in all lakes, and with strong effect sizes. Its inclusion in other focused tests reduced the effect size of the modular signal, particularly for the defensive trait suite, indicating shared variation between defensive traits and some components of shape (Table S14).

Responses to the environment
The surveyed lakes are mostly of small to moderate size (min. surface area = 2.7 ha, Ruth Lakealthough Ruth was excluded from PCAs due to lack of Ca data; max. surface area = 135.7 ha, Finger Lake). However, these lakes exhibit a wide range of environmental characteristics.
Lakes ranged from 2.1 (Corcoran) -24.4 m (Wik) in maximum depth. All lakes had low to very low biomass of primary producers (max. Chl-a = 3.30 µgL -1 , Finger Lake; min. Chl-a = 0.42 µgL -1 , Wik Lake). The pH in the lakes ranged from circumneutral (pH = 6.6, Ruth Lake) to Figure 6. PCA biplots of lakes and environmental variables. Panel A shows physico-chemical environmental variables. Panel B shows a PCA biplot additionally including Daphnia abundance, Gammarid abundance and percent littoral area, all foraging-related environmental variables. Ruth Lake and Echo Lake are excluded from both panels because some environmental data were unavailable.

A B
alkaline (pH = 8.8, Corcoran Lake), and Ca and conductivity were both very low in some lakes.
The two lakes with the lowest Chl-a, Wik and G Lakes, also had the lowest Ca (1.3 and 0.6 mgL -1 , respectively) and cond. (14.9 and 11.3 µScm -1 , respectively). Jean Lake had the highest Ca  Two-block partial least squares regressions between trait suites and environmental variables did not reveal any significant correlations after p-values were adjusted for multiple comparisons (Figure 7; Table S17). However, physico-chemical environmental variables alone were correlated with stickleback defensive and swimming trait suites at levels only slightly higher than the p = 0.05 threshold. While the trophic trait suite was not strongly correlated with physico-chemical characteristics of the lakes, this association was strengthened substantially by the inclusion of environmental variables related to foraging (Daphnia sp. abundance and gammarid density, and % littoral area), even though this reduced the number of lakes for which data was available to only eight lakes, and the correlation remained statistically insignificant.
Axis loadings and singular values for all PLS tests are shown in Tables S18-S21. When the angles of the primary environmental vectors were calculated between 2B-PLS tests against different trait suites, angles involving the stickleback defensive trait suites were consistently largest. The defensive trait suite angles, especially the angles between the defensive trait suite environmental axis and the trophic trait suite environmental axis, indicate that the major axis of divergence in defensive traits among lakes occurs along an environmental gradient that is divergent from those of the other traits (Figure 7; Table 3). These angles also widen greatly with the inclusion of foraging-related environmental data, further emphasizing the distinction between the axis of divergence in defensive traits and the more congruent trophic and swimming traits.

Discussion
Even on short timescales, stickleback populations can adapt in response to selection on several trait types simultaneously (Marques et al. 2018). However, because traits often covary, even when their functional relationships are not obvious or are indirect (Barrett et al. 2008;Bjaerke et al. 2010), it has been difficult to establish whether populations diverge along a lowdimensional space following a genetic line of least resistance (Schluter 1996), or diverge in a higher-dimensional space as a consequence in response to multifarious selection on functionally disparate traits (Nosil et al. 2009).
We show that divergence is not a unidimensional process, but instead takes place along several axes both within and across trait suites under selection from distinct sources. This permits multiple possible evolutionary trajectories through phenotype space depending on the selection pressures experienced. While some within-lake patterns of integration between suites of traits are consistent across lakes, others vary considerably in both strength and direction between populations. Consistency and extent of trait suite modularity also varies by trait suite. The major axes of variation in the swimming and trophic trait suites also correspond to PLS axes of similar environmental variation. However, these diverged from the environmental axis along which defensive traits varied, indicating potential tradeoffs in which the optimum phenotypes of different trait suites diverge.
The most multidimensional patterns of divergence existed in body shape data and in the trophic trait suite. For the defensive trait suite, there was low dimensionality of divergence, but this was partially due to the finding that divergence along LD1, dominated by pelvic spine length, was so great that it swamped other dimensions in the calculation, rather than the lack of divergence along other axes. Importantly, while primary axes of divergence were strongly correlated between trait suites, trait suites were less strongly correlated with one another when they were considered as multivariate trait matrices. These patterns illustrate that lowermagnitude components of divergence may be more susceptible to selective pressures not shared between trait suites. This conclusion is also supported by the persistent, but relatively smalleffect overall signal of modularity between trait suites.

Patterns of integration between trait suites often vary in strength and direction between
populations, demonstrating that important axes of divergence not only exist in the trait values themselves, but also in the ways that trait suites are related to one another. This could potentially influence how trait suites respond to the environment via both adaptation and plasticity. In addition to our conclusions about the dimensionality of divergence within and between trait suites, we established that the defensive trait suite was the most independent of the trait suites we considered. This was evident from the inconsistent orientation of its axes of covariation with the other trait suites, and also from its high levels of modularity in all lakes, when considered independently of shape. 2B-PLS tests of trait suites against environmental variables further showed that the environmental gradient along which defensive traits varied diverged in its orientation from the environmental gradients along which the other trait suites varied. Trophic and swimming traits varied much less independently than defensive traits, although some lakes did exhibit modularity in trophic phenotypes. The axes along which they were correlated with one another were also very similar among lakes.

General Conclusions
Three notes of caution regarding dimensionality of divergence should accompany our findings. 1) The range of habitats that we have included, while ecologically diverse for lakes with stickleback populations, represent only a small subset of the range of habitats in which stickleback are found, and consequently, probably underestimate the magnitude, and likely the dimensionality, of divergence among stickleback populations more broadly. There are also anadromous marine stickleback populations, as well as stream populations (Lavin and McPhail 1993;Bell and Foster 1994;Seebacher et al. 2016). Inclusion of either or both population types would likely have substantially increased the dimensionality calculated for the defensive trait suite, and the shape and swimming trait suites, respectively. In the former case, this is because threespine stickleback have three major lateral plate morphs: low, partial, and complete, and while low-and partial-plated morphs are common in freshwater, marine populations are typically made up almost entirely of completely-plated morph (Hagen 1967;Bell and Foster 1994;Barrett et al. 2008). Of the stickleback represented in this study, there was only one single partial-plated individual (from Spirit Lake), and the remainder were low-plated. This means that the fish we sampled only represent a small band of the possible variation in this trait, and divergence in defensive traits was determined primarily by spine length in our data. In the latter case, populations of stickleback inhabiting streams typically have trophic adaptations suited to benthivory, as benthic invertebrates are more reliably available as prey items in most streams than zooplankton (Berner et al. 2008;Berner et al. 2009). They also typically have a deeper body shape that facilitates navigating the spatially heterogeneous flow regimes present in many streams (Hendry and Taylor 2004;Izen et al. 2016). The specificity of our study to lacustrine populations likely does not adequately represent the dimensionality of trait divergence across the entire species or the diversity of integration structures among trait suites.
2) Another important consideration is that our results seem to indicate that both effective dimensionality and variability in the orientation of integration structure PLS axes generally increase with the number of traits under consideration, and this observation was supported by our rarefaction-like analysis in the appendices. This may be intuitiveafter all, our swimming trait suite, which is composed of only three traits, cannot possibly have more than three effective dimensions. Neither can the number of effective dimensions of divergence exceed one fewer than the number of populations under study. Nevertheless, the intrinsic dimensionality (Del Giudice 2020) of divergence seems to not often be considered when authors assess whether an adaptive radiation is constrained to a line of least resistance (Schluter 1996). It is therefore important that researchers, when testing and making claims about the dimensionality of adaptation or divergence, measure broad ranges of traits of differing functions in populations of diverse habitat types, unless the intention is to test a narrow question relating to patterns of adaptation in a particular suite of traits. A rarefaction-like procedure similar to the one we have implemented in the appendices could help fewer researchers avoid the hazards involved in calculating effective dimensionality.
3) Finally, we advise other researchers to be conscious not just of the dimensionality of divergence in their systems, but also of the magnitude of divergence. The measures of effective dimensionality we have employed represent relative strength of the vectors of divergence to one another, but do not describe the extent to which they have diverged. In our case, although the populations we sampled are visibly divergent without measurement of any phenotypes, the large magnitude of the difference between the populations with pelvic spines and those without pelvic spines dwarfed other dimensions of divergence. Spine length reduced effective dimensionality of traits, and especially of the defensive trait suite, even when differences in other traits among populations were substantial, and this effect can be seen in Figure B1. Consideration of the magnitude is also necessary for measurements of divergence to be useful as adaptive responses to the environment (Stuart et al. 2017;Runquist et al. 2020), or when assessing the effect of divergence on an environment.
Our findings suggest that divergence among populations is not typically constrained by the axis of greatest covariation, even within groups of traits that might be expected to respond to similar selection pressures. This conclusion is consistent with the known examples of highly divergent conspecific sympatric morphs even within a single habitatespecially among salmonids (Skúlason et al. 1989;Jónsson 2000;Chaverie et al. 2016;Muir et al. 2016;Arostegui and Quinn 2019;Chaverie et al. 2020), and is incompatible with a simple singly-peaked orridged fitness landscape. Although evidence from recent post-glacial radiationslike some sticklebacks and the Salvelinus speciesare particularly compelling because of the relatively short period of time available for the reorientation of gmax in a low-dimensional model of divergence, other radiations also provide support for multi-axis divergence that is inconsistent with low-dimensional adaptation constrained by genetic covariances (Velasco and Herrel 2007;Hulsey et al. 2019;Levin et al. 2021). Further investigations into divergence as it progresses could further elucidate this process. This is especially true of the role of plasticity, which Noble et al. (2019) have shown is often aligned with gmax, which would be expected to facilitate lowdimensional divergence, in spite of our results.
The implications of these findings are important to consider when attempting to assess the repeatability and predictability of evolution in nature, one of the central themes of evolutionary biology. The possibility of adaptation in a high-dimensional space certainly makes predicting the evolutionary effect of an environmental condition in any particular instance difficult, even with high-quality data about multiple aspects of the environment and when the relatedness among populations is known (Stuart et al. 2017;Yong et al. 2020). However, the potential to adapt along a number of axes may increase resilience to environmental changes, allowing a population to adapt to the change in the context of its specific habitat conditions (Reisch et al. 2015), rather constraining it to a single axis that may be perpendicular to the direction of selection. Because the dimensionality of potential divergence is likely very taxonspecific, further investigations into the genomic mechanisms and the role of plasticity through time in diverging populations will be necessary to develop useful predictive models for the evolution of species that are vulnerable, ecologically important, or important to human health. Table A1. Tests in the Methods, and the Panels of Figure 1

to which they Correspond Question
Test Figure 1 panel (

The Use of Linear Discriminant Analysis for Calculation of Effective Dimensionality
We use LDA, which places an emphasis on the differences between groups relative the variance within groups, rather than Principal Components Analysis (PCA) to identify axes of divergence because the primary subject of our analyses is the dimensionality of divergence between populations, not simply differences between group means. Typically dimensionality of divergence is studied using eigenvectors calculated from population means (as in Schluter 1996), but because of differing magnitudes and orientations of trait variances within populations, this does not capture the axes along which populations are actually most divergent from one another.
As such, the use of eigenvectors of population means represents the process of phenotypic divergencewhich can be constrained by the evolvability of traits, lineage history, and the ability of variation to persist in populations due to population size and balancing or disruptive selectionas equally possible in all directions. Effectively, it is a mathematical convenience that is agnostic towards which constraints on phenotypic change are most important (which may differ depending on study question).

Comparison of Effective Dimensionalities of Trait Suites with Differing Numbers of Traits
Effective dimensionality of data matrices depends in part on the number of traits included in the data. For instance, the data we include in our swimming trait suite includes only three traits, and it therefore not possibly have an effective trait dimensionality higher than three. As a consequence, it was necessary to develop a method that would enable comparing trait dimensionality between trait suites at different numbers of traits. We accomplished this by implementing a rarefaction-like procedure in which traits were randomly drawn without replacement from the trait matrix for each trait suite. After each draw, an LDA was performed with the sampled traits, and dimensionality of divergence calculated using the Cangelosi and Goriely / Roy and Vetterli index presented in Table 1

Diet-Related Environmental Data
To estimate Daphnia abundance in the pelagic zone of each lake, crustacean zooplankton were collected by whole water column vertical tows with a 35-cm diameter Wisconsin net from 0.5 m off the bottom to the surface. The zooplankton were sampled from four sampling stations along an open-water transect across each lake. The crustacean zooplankton were anaesthetized with bromoseltzer, preserved with 95% ethanol, and they were brought back to the laboratory at the UQAM for identification. These four samples were subsequently pooled to a single sample per lake for identification and enumeration in the laboratory. The crustacean zooplankton were identified to species level with a high-resolution dissecting microscope (6.3-126x; SZ2-IL-ST; Olympus SZ) and species were enumerated for the number of individuals per L. Crustacean zooplankton were counted using a protocol that targeted mature individuals that could be identified unambiguously to species, as well as to detect rare species (Mack et al. 2012).
Subsamples (10 mL) were taken from a standardized 50 mL sample volume, until at least 200 individuals of the most abundant species were counted (excluding copepodids and nauplii), or a total of 1000 individuals were counted (excluding copepodids and nauplii), or the total sample was counted. Taxonomic keys for Daphnia sp. included Brooks (1957), Taylor and Hebert (1993), Witty (2004), Thorp and Covich (2010), and Haney et al. (2013).
To estimate gammarid density in the littoral zone of each lake, nearshore littoral macroinvertebrate communities were collected with a D-frame kick net at four sampling stations located around the lake in the 2018 survey, and eight sampling stations in the 2019 survey.
Sampling was performed by the sweep method ("Kick & Sweep") with a 500 μm "D-net" as recommended by the Ontario Benthos Biomonitoring Network (Jones et al. 2007) on a surface of approximatively 2 m 2 . Samples were concentrated with a 500-m sieve, preserved in 95% ethanol, and they were brought back to the laboratory at UQAM for identification. Identification was done up to the Family level following taxonomic keys (Merritt et al. 2008;Moisan 2010).
All macroinvertebrate samples were identified using SZX10 stereo microscopes (Olympus) with varying magnification (x6.3 -x10). For each sample, 100mL sub-sample was taken and counted until the 100th individual was reached. If the 100th individual was not reached within the initial 100 mL, an additional 100 mL was counted. When the 100 th individual was reached, the remaining part of the sub-sample was counted, and the total sub-sampled volume calculated. A ratio was then calculated between the total sub-sampled volume and the sample total volume to estimate the taxon-specific abundance of macroinvertebrates. The density of gammarid individuals per m 2 was averaged for the 4 samples collected per lake to obtain total average gammarid density per lake in the littoral zone.

Percent Littoral Area
To calculate the proportion littoral area for each lake, we used bathymetric maps that were available through the Alaska Department of Fish & Game. For each of the maps' contours, we calculated the area by measuring the area of the map image bounded by the contour in ImageJ (Rasband 1997(Rasband -2018. This was then converted to area using a scale indicated on the map, a total surface area indicated on the map, or, when no indication of scale was indicated on the map, the total surface area from previous surveys used in the environmental data. When maximum depth was not indicated on a map, the maximum taken during field surveys was used, unless there was a contour on the map below this depth. For Tern Lake, the only bathymetric map available used transects, rather than contours, to indicate depth. In this case, transect intervals were assumed to be proportional to area, although this likely resulted in an overestimate of the area of the deeper portions of the lake. The total area of Tern Lake as indicated on the documents accompanying the bathymetric map was also much smaller than the area used in the Because bathymetric maps from the different lakes had depth contours marked at inconsistent intervals, it was necessary to estimate depths at common intervals across all lakes. Therefore, after collecting depth profile data from bathymetric maps, we computed a monotonic spline for each lake using surface area as a response variable to depth. The splines used the method of Hyman (1983) as implemented in R's stats package (R Core Team 2019), to avoid the calculation of increases in area with increasing depth. The littoral area was defined as the area of the lake having a depth of three meters or less, and the percent littoral area was calculated as 1-AreaDepth>3m / AreaTotal where the area with a depth > 3 m was interpolated from the spline.

Comparison of Effective Dimensionalities of Trait Suites with Differing Numbers of Traits
The permutation procedure implemented to compare effective dimensionalities of differing trait groups showed that at low numbers of traits, the effective dimensionality of swimming and trophic trait suites is slightly lower, but similar, to that of shape ( Figure B1).
Divergence in the defensive trait suite is consistently much lower in effective dimensionality, as is divergence in the full non-shape trait data. However, this is primarily due to the inclusion of pelvic spine length, which results in the clear decline in divergence dimensionality in a permutation once it is included as a trait, an effect visible in the Defensive and Full data in Figure B1. Prior to inclusion of this trait in a permutation, the effective dimensionalities for both the defensive and full non-shape data are comparable to those of shape, trophic, and swimming trait suites. Figure B1. Effective dimensionality of divergence for each trait suite, as well as the full non-shape phenotype data, with the inclusion of different numbers of traits. Thin lines represent 50 permutations of randomly sampling traits until all traits are included in the dimensionality calculation of LDA vectors. The left panel includes complete data from the permutation procedures, and the right panel shows only the first five traits in the permutations along with loess regressions showing the mean effective dimensionality for each trait suite at a given number of traits.

Diet-Related Environmental Data
Several patterns emerged from diet-related environmental data. For the Daphnia abundance, lakes clustered into two groups: Watson, Wik, and Spirit Lakes had over 5 Daphnia individuals L -1 , and all other lakes, including Corcoran and Ruth (which are not plotted due to the absence of littoral area data), had less than 0.15 Daphnia L -1 ( Figure B2). Most lakes had gammarid abundance less than 10.5 individuals per m 2 , with Long Lake having 10.5 individuals per m 2 , and Wik Lake having the fewest, with 1.0 individuals per m 2 . Tern Lake, however, had the highest abundance of gammarids, with 22.8 individuals per m 2 . Of the lakes with the highest percent littoral area -Tern (75.0%), Watson (75.6%), and Walby (82.1%) -Tern and Watson fish correspondingly appeared at the "benthic" extreme on the primary trophic trait LDA axis, although Walby Lake was more moderate in trophic morphology.
Interspecific competition likely played a role in the unusual prey abundance data for Watson Lake, which has threespine stickleback morphologies that we would expect of a benthivorous population, despite the abundance of Daphnia and paucity of gammarids. Of all the lakes sampled, Watson Lake was the only lake with a large ninespine stickleback (Pungitius pungitius) population. Wik Lake was the only other lake in which any ninespine stickleback were caught in 2018, but this was only a single individual. In a subsequent 2019 survey of Wik Lake, no ninespine stickleback were captured. State records also indicate that the two lakes immediately downstream of Watson, which are connected by short lengths of stream easily passable by fish, also have abundant sizeable populations of rainbow trout (Oncorhynchus mykiss) and longnose sucker (Catostomus catostomus) (ADF&G 2015). These species all likely compete with threespine stickleback for benthic prey, such as gammarids, and although we did not have sufficient data to formally test this the eco-evolutionary effects of competition here, this is a potentially fruitful avenue of future work.

Correlations of Environmental Variables with Phenotypes
Spearman's correlations revealed that a large number of traits are associated with Calcium, pH, and conductivity, which are highly correlated environmental characteristics in these lakes ( Figure B3). These traits include plate number and spine length, which are expected to be associated with are positively associated with Ca, pH, and conductivity, due to the dependence of these skeletal features on available calcium (Bell et al. 1993;Klepaker et al. 2013; Figure B2. Density of Daphnia and Gammaridsimportant crustacean prey items for zooplanktivorous and benthivorous stickleback, respectively. Daphnia abundance is shown as individuals per L, and Gammarid abundance is shown as individuals per m 2 . Color scale indicates the proportion of each lakes surface area shallower than 3 m. Klepaker et al. 2016 with the cost of high-platedness on fast-start swimming (Bergstrom 2002).
In terms of the trait suite discriminant axes, the primary axis of defensive trait divergence is negatively associated with nutrient and Ca ion availability, as well as benthic prey abundance.
The primary axes of trophic and swimming trait divergence are both strongly associated with lake depth and DOC concentrations, with the more "benthic" ends of these axes associated with A B Figure B3. Spearman's correlations between environmental variables and stickleback traits (A) and trait suite discriminate axes (B). Color scale indicates Spearman's .
shallow, low DOC lakes. They are associated with different prey types, however, with more "benthic" swimming values negatively associated with Daphnia abundance, and more "benthic" trophic values counterintuitively negatively associated with gammarid abundance. Despite the previously noted correlations in the LD axes, therefore, their environmental associations make clear some incongruities in the environmental correlates of divergence among trait suites.