## Abstract

Phenotypic (co)variation is a prerequisite for evolutionary change, and understanding how (co)variation evolves is of crucial importance to the biological sciences. Theoretical models predict that under directional selection, phenotypic (co)variation should evolve in step with the underlying adaptive landscape, increasing the degree of correlation among co-selected traits as well as the amount of genetic variance in the direction of selection. Whether either of these outcomes occurs in natural populations is an open question and thus an important gap in evolutionary theory. Here, we documented changes in the phenotypic (co)variation structure in two separate natural populations in each of two chipmunk species (*Tamias alpinus* and *T. speciosus*) undergoing directional selection. In populations where selection was strongest (those of *T. alpinus*), we observed changes, at least for one population, in phenotypic (co)variation that matched theoretical expectations, namely an increase of both phenotypic integration and (co)variance in the direction of selection and a re-alignment of the major axis of variation with the selection gradient.

## 1. Introduction

In order to persist over evolutionary time, species must have the ability to respond in the direction of selection. As organisms are formed by a combination of multiple traits organized into a coherent whole (a multidimensional system), understanding the interaction between the available phenotypic (co)variation and selection is crucial to understand species' responses to selection [1] and, consequently, species' persistence over time. For instance, if a species lacks phenotypic (co)variation in a certain direction it can quickly become extinct when directional selection operates along that trajectory [2]. How the available phenotypic (co)variation shapes species evolution is therefore an important area of research in evolutionary biology, and is relatively well appreciated theoretically, empirically and computationally (e.g. [1,3,4]). On the other hand, our understanding of how directional selection shapes the evolution of available phenotypic (co)variation, although equally important, is only just beginning. Most studies addressing this issue have been simulations derived from theory [5,6], although some experimental evidence has emerged [6–9].

In traditional evolutionary thinking, directional selection is thought to deplete genetic variation, leading to a decrease in phenotypic (co)variation [2,10]. Therefore, establishing plausible mechanisms that account for the widespread phenotypic (co)variation observed in nature became a priority in evolutionary biology [10]. Despite the inherent difficulties in pursuing answers to this issue [5], recent advances in theoretical and computational studies have provided some benchmarks for how available phenotypic (co)variation is expected to evolve under directional selection [11–15]. These models are based on the evolution of the genotype–phenotype map (GP-map), which describes how genetic (co)variation is translated into phenotypic (co)variation. If different genotypes differ in genetic (co)variation among traits (i.e. the degree of pleiotropy and epistasis among traits), this leads to the possibility that genetic (co)variation among traits can evolve in response to selection [13,14]. In these models, selection might actually drive the evolution of the mutational machinery to align available phenotypic (co)variation with selection [16]. Consequently, the impact of directional selection on the GP-map organization is thought to be substantial and can occur at a rapid pace [12,13]. Such models allow us to make predictions about how change occurs in the available phenotypic (co)variation that are easily testable by empirical data. First, we would expect an increase in genetic variation in the direction of selection (figure 1). Second, the degree to which traits are correlated is expected to change; specifically, we expect increased correlation among traits that are being co-selected [6]. Lastly, we expect a re-orientation of the pattern of available phenotypic (co)variation to match the expected direction of selection [12] (figure 1).

Validation of these predictions in multidimensional traits from natural populations is essential to further develop our understanding of evolution itself and of how species adapt to new selective pressures. It is especially relevant in a world where most natural environments are under some kind of stress due to anthropogenic pressures, either direct (e.g. land use change) or indirect (e.g. climate change). However, three main problems hinder the gathering of empirical data to assess these questions in natural populations. First, sampled populations should be separated by many generations, given that such changes are not expected in a short period of time. Second, populations must be well sampled in order to properly estimate statistical parameters. Third, the populations in question should have experienced directional selection. Here, we used a unique sample set that matches the first and second requirements, and which thus permits us to test the third requirement.

Our sample is composed of chipmunk specimens of two species, *Tamias alpinus* and *T. speciosus*, collected almost a century apart. These come from two independent transects separated by approximately 180 km along the Sierra Nevada of California. The first of these is in the central part of the Sierras, in Yosemite National Park; the other is in the Southern Sierras, within or just east of Sequoia-Kings Canyon National Parks. These species are phylogenetically close, but have responded in strikingly different manner to a century of climatic and associated environmental changes. While the alpine chipmunk, *T. alpinus,* has shifted its elevational distribution upwards and decreased its genetic diversity, the lodgepole chipmunk, *T. speciosus*, has changed neither its elevational distribution nor genetic diversity [17–19].

In this study, we tested for both species whether evolution between the historical and modern periods was driven by directional selection, and because we found evidence favouring a directional selection scenario, we tested how the selective regime changed the available phenotypic (co)variation over time. To this end, we analysed 35 skull traits in those populations using a quantitative genetics framework. We then compared the phenotypic (co)variation matrices between periods (historical and modern) for each species and transect in order to assess whether the specific selective regimes in each population had an impact on the overlying phenotypic (co)variation patterns.

## 2. Material and methods

### (a) Samples

One of us (A.P.A.A.) recorded three-dimensional coordinates for 27 landmarks on 197 adult skulls of *T. alpinus* (Yosemite: 51 historical/38 modern; Southern Sierras: 75 historical/33 modern) and 481 adult skulls of *T. speciosus* (Yosemite: 77 historical/221 modern; Southern Sierras: 83 historical/100 modern). Adult specimens were defined by a fully erupted permanent premolar 4 and a completely fused basisphenoid–basioccipital suture. Based on these landmarks we calculated 35 linear distances, which were used in the subsequent analyses (figure 2). We chose to use linear measurements instead of landmark coordinates because we are primarily interested in differences in covariance structure. Proper estimation of the covariance structure in a geometric morphometrics approach would require much larger sample sizes than were available (at least three times the number of landmarks) [20].

We used historical samples collected as part of a California-wide survey of terrestrial vertebrates conducted by Joseph Grinnell and colleagues from 1911 to 1915. Modern specimens were collected as part of the Grinnell Resurvey Project, an intensive resampling of Grinnell's historic sites that occurred from 2003 to 2013, and were collected under the auspices of the US National Park Service, specifically Yosemite and Sequoia-Kings Canyon National Parks, and under permit from the California Department of Fish and Wildlife, the Animal Care and Use Committee of the University of California at Berkeley, and followed the guidelines of the use of wild mammals in research developed by the American Society of Mammalogists [21]. All specimens are deposited in the Museum of Vertebrate Zoology at the University of California, Berkeley. Samples come from two independent transects, one through Yosemite National Park in the Sierra Nevada in central California and the other from the Southern Sierras, within and just east of Sequoia National Park approximately 180 km to the south. It has been previously shown [22] that the *T. alpinus* population from Yosemite has changed its skull morphology in a pattern compatible with directional selection when tested in a univariate context, with most changes being concentrated in the facial region for both transects. Also, while the Yosemite transect population has increased in size over the past century, the Southern Sierras population has decreased in size [22]. By contrast, comparable samples of *T. speciosus* had fewer traits that changed in a pattern compatible with directional selection in either transect, although in Yosemite most changes were also concentrated in the facial region.

### (b) Matrix similarity

Here, we used the phenotypic matrix (**P**-matrix) as a proxy for its genetic counterpart, which is the evolutionary relevant parameter. Our decision to use **P**-matrices as a substitute for **G**-matrices is based on considerable evidence supporting their interchangeability, at least for morphological characters, especially in mammals [23–25]. Furthermore, to guarantee that our **P**-matrices were similar to other estimates of **G**, we compared our **P**-matrices with a *Calomys expulsus* **G**-matrix [26]. Also, we evaluated the overall similarity between the covariance matrices of all populations (historical and modern) using both the random skewers and the Krzanowski methods [27,28]. The random skewers method is based on the multivariate response to selection equation. This method consists of applying a set of random selection vectors (normalized to a length of one) to a pair of matrices and comparing the response vectors between the two matrices using a vector correlation procedure. The mean vector correlation between the matrices' evolutionary-response vectors is a measure of similarity. The more similar the response to selection, the higher the similarity between the matrices [28]. The Krzanowski method consists of comparing a pair of matrices by calculating the angles between the best-matched pairs of orthogonal axes (principal components) [27,29] (see the electronic supplementary material for further discussion). As matrices are estimated with error, we corrected matrix similarity values using matrix repeatability estimates [28] (electronic supplementary material, table S1).

### (c) Directional selection versus genetic drift of skull traits

Before exploring whether directional selection had an impact on available phenotypic (co)variation, we tested whether natural selection was responsible for multivariate evolution between the historical and modern periods. To do so, we used the theoretical expectation that the amount of divergence expected for a population evolving under drift, in a multivariate framework, corresponds to a divergence matrix **D**. An estimation of this **D**-matrix can be made by **D** = **G**t/Ne, where **G** represents the **G**-matrix of the ancestral population, t is the divergence time in generations, and Ne is the effective population size [1,30]. We estimated the **D**-matrix using the historical **P**-matrix as surrogate for the ancestral **G**-matrix (but see the electronic supplementary material for estimates using different matrices, table S5), Ne derived from upper and lower estimates for *T. alpinus* [31], and a generation time of 1 year [32]. Having estimated the **D**-matrix, we used it as the *Σ* parameter in a multivariate normal distribution, from which we randomly sampled 1000 vectors. Each of these vectors represents potential trait values for a population evolving under drift. The subtraction of each of those 1000 random vectors from the ancestral populations' trait means gave us estimations of morphological change expected under drift (Δ_{z} expected in a drift scenario). Finally, we compared the 95% probability distribution of the magnitude of morphological change expected under drift (which was estimated as the norm of the randomly sampled Δ_{z} vectors) with the empirical norm of Δ* _{z}* observed (i.e. estimation based on the historical and modern populations data, taking into account mean standard errors). We then concluded that directional selection was the main process responsible for observed divergence if the range of estimates from the magnitude of morphological divergence fell outside the 95% distribution expected under drift.

### (d) Selection gradient estimate

In order to estimate the observed net selection gradient (*β*), which is the directional selection responsible for the observed morphological changes, we used the multivariate selection response equation [1] where Δ* _{z}* is the response to selection estimated as difference between modern and historical observed trait means;

**G**is the inverse of the

^{−1}**G**-matrix (substituted by the historical

**P**-matrix). Owing to the fact that inverted matrices are dominated by small eigenvalues usually estimated with a lot of noise, we controlled our

**P**-matrices for noise using an eigenvalue extension method [33]. The calculated normalized

*β*was then used as a benchmark to estimate changes in the pattern of (co)variation, as described below. Moreover, to compare the strength of selection between transects and species, we estimated the norm of the selection gradient standardized by trait means [34].

### (e) Effects of directional selection on skull **P**-matrices

As our main purpose was to compare changes in the pattern of (co)variation between time periods, we had to understand the possible distribution of those parameters in each species/area sample. To do so, we estimated **P**-matrices separately for the historical and modern samples, estimating a distribution of 1000 **P**-matrices by Monte Carlo resampling for each period, transect and species [35]. We then used these 1000 estimated matrices in subsequent analyses and considered a result significant whenever the observed modern **P**-matrix value fell outside the 95% distribution of the historical estimates [36]. To circumvent potential biases in our comparisons due to differences in sample sizes, the random populations sampled in the Monte Carlo procedure were set to the lowest sample size between the time periods compared.

To quantify the impact of natural selection on the **P**-matrices, we compared historical and modern **P**-matrices in relation to three different matrix features: (i) size, which can be described as the total amount of variation in the matrix or in a certain direction [34,37]; (ii) shape, which provides an indication of eccentricity, or how tight the correlations between traits are; and (iii) orientation in relation to the selection gradient.

In order to determine if the total amount of variance had changed from historical to modern samples, we estimated the trace of each covariance matrix from the Monte Carlo distribution. In addition, we determined if the amount of variance had changed in the direction of selection or in directions uncorrelated with the observed selection gradient. The amount of variance in any given direction was calculated as the evolvability in that direction [34], a metric which captures the ability of a population to evolve in the direction of a specified selection gradient. Evolvability can be measured as
where e(*β*) is the evolvability in the direction of a given selection gradient (*β*) and **G** is the **G**-matrix. To compare the effects of directional selection, we estimated the evolvability in the direction of the normalized observed selection gradient (*β*_{obs}) for the distribution of **P**-matrices from the historical and modern periods. Moreover, we generated 1000 random selection gradient vectors uncorrelated with *β*_{obs}. In order to obtain those sets of uncorrelated vectors, we first sampled 1000 normalized vectors from a normal distribution and applied the formula
where *β*_{r} is the random selection gradient sampled from a normal distribution and *β _{i}* is the uncorrelated resulting vector. We later normalized the random vectors to 1 and compared the evolvability potential of both the historical and modern

**P**-matrices in those directions. Evolvability estimates were divided by the geometric means of all traits, thus accounting for scale differences between populations [34].

We also compared changes in the overall magnitude of integration between periods by estimating the coefficient of determination (*r*^{2}) of the correlation matrices [38]. This coefficient is simply the average of the squared correlation coefficients and measures how tightly the correlations among traits are, with a higher estimate indicating higher correlations.

Lastly, we compared changes in the orientation of the **P**-matrix distributions between periods by estimating the correlation between the observed selection gradient and the first principal component of each of 1000 matrices from historical and modern times for both species and transects. All statistical analyses were done in the R environment for statistical computing [39] using the EvolQG package [35].

## 3. Results

### (a) Matrix comparisons

Phenotypic matrices (**P**-matrices) were used as proxies for their underlying genetic matrices (**G**-matrices), which are the evolutionarily relevant parameter. Consequently, prior to any analyses we had to ensure that our estimated **P**-matrices were good estimates of the underlying **G**-matrices. To this end, we compared our **P**-matrices with a **G**-matrix derived from the distant related rodent *Calomys expulsus* [26] using both the random skewers and Krzanowski methods [8,28,29]*.* These comparisons showed a high degree of similarity between the **P**-matrices and the **G**-matrix, with comparisons from the random skewers method ranging from 0.61 to 0.81 (all significant, indicating that they are more similar than would be expected by chance) and those from the Krzanowski method ranging from 0.66 to 0.72 (electronic supplementary material, table S2). Although there is no consensus among researchers as to how similar a matrix has to be to be considered similar [23,40,41], these fairly high values indicated that our **P**-matrices were reasonable estimates of their genetic counterparts [24,41]. Moreover, comparisons of the **P**-matrices among populations showed that all were structurally similar, with random skewers estimates ranging from 0.81 to 0.95 and Krzanowski estimates ranging from 0.81 to 0.87 (electronic supplementary material, table S3).

### (b) Directional selection versus genetic drift of skull traits

To determine whether the amount of divergence observed between historical and modern periods for each population was explained by genetic drift or by natural selection, we simulated the amount of morphological divergence expected by drift and compared it with the empirically measured magnitude of morphological change. For any effective population size adopted, the magnitude of empirically measured morphological change in the Southern Sierras was 5.7 and 10.95 times higher than that expected by drift for *T. speciosus* and *T. alpinus*, respectively. For the Yosemite transect, the same pattern was observed; the magnitude of empirically observed morphological change was 7.5 and 19.6 times higher for *T. speciosus* and *T. alpinus*, respectively, than that expected by drift. Therefore, data for both species support directional selection as the primary mode underlying the observed temporal changes (table 1).

Next, we estimated the standardized magnitude of morphological change (||Δ*zµ*||) and selection gradient (||ß*µ*||), estimated as the norm of each vector, for both species to gauge the strength of selection. For *T. speciosus* from the Yosemite transect, we obtained ||Δ*zµ*|| = 0.064 and ||ß*µ*|| = 15.457; comparable numbers for *T. alpinus* were ||Δ*zµ*|| = 0.178 and ||ß*µ*|| = 39.388, which were 2.8 and 2.6 times higher, respectively, than those of *T. speciosus*. In the Southern Sierras, we obtained ||Δ*zµ*|| = 0.063 and ||ß*µ*|| = 20.279 for *T. speciosus*, while the results for *T. alpinus* were ||Δ*zµ*|| = 0.099 and ||ß*µ*|| = 27.062, values that were 1.57 and 1.33 times higher, respectively, than those of *T. speciosus*. In both cases, directional selection was stronger in *T. alpinus* than in *T. speciosus*.

### (c) Effect of directional selection on morphological **P**-matrices

To determine the effect of directional selection on the **P**-matrices, we investigated the following changes between the historical and modern **P**-matrices: (i) the total amount of variation estimated by the matrices' traces; (ii) the amount of variation associated with the direction of selection and with directions uncorrelated with selection, calculated as evolvability divided by the geometric mean of the traits; (iii) the overall magnitude of correlation among traits estimated as the coefficient of determination of the correlation matrices (*r*^{2}); and (iv) the orientation of the axis of greatest variation in relation to the selection gradient estimated by the correlations between PC1 and the selection gradients.

For *T. alpinus* from the Yosemite transect, the historical and modern matrix traces were 2.79 ± 0.22 s.d. and 3.33 ± 0.34 s.d., respectively, while for *T. speciosus* the corresponding values were 4.34 ± 0.28 s.d. and 4.12 ± 0.23 s.d., respectively. For the Southern Sierras transect, we obtained historical and modern matrix traces for *T. alpinus* of 2.63 ± 0.24 s.d. and 2.60 ± 0.25 s.d., respectively, while the corresponding values for *T. speciosus* were 4.68 ± 0.30 s.d. and 4.47 ± 0.31 s.d, respectively. Thus, the total amount of variation in each species-population comparison did not change temporally. However, the amount of variation in the direction of selection did increase for *T. alpinus* in both transects (figure 3*b*; Yosemite observed modern estimate = 0.033, 95% historical distribution = 0.0128–0.031; Southern Sierras observed modern estimate = 0.018, historical distribution = 0.006–0.015) but not for *T. speciosus* (figure 3*d*; Yosemite observed modern estimate = 0.024, 95% historical distribution = 0.015–0.029; Southern Sierras observed modern estimate = 0.011, historical distribution = 0.009–0.017). Importantly, there was no change in the amount of variation in either species or transect in directions uncorrelated with the selection gradient (figure 3*a,c*).

The overall magnitude of integration increased over time for *T. alpinus* in Yosemite, as the observed *r*^{2} index for the modern population (0.108) did not overlap with the historical distribution (0.072–0.106). For *T. alpinus* from the Southern Sierras, however, the overall magnitude of integration remained unaltered across time (observed modern *r*^{2} index = 0.094, 95% historical distribution = 0.079–0.140). We also observed idiosyncratic changes in *T. speciosus*, as the Yosemite population decreased its overall magnitude of integration (the observed modern value of 0.073 does not overlap the historical 95% distribution of 0.078–0.131), while the Southern Sierras population remained unaltered between the historical and modern periods (observed modern value = 0.118, 95% historical distribution = 0.085–0.133; figure 4).

Lastly, the orientation of the axis of greatest variation in relation to the selection gradient, estimated by the correlations between PC1 and the selection gradients, did not change for most populations (figure 5). The only population where we observed an increase in the correlation between PC1 and the selection gradient was *T. alpinus* from the Southern Sierras, where the modern observed correlation of 0.31 was larger than the historical 95% distribution of 0.027–0.23.

## 4. Discussion

How species adapt to their environment is a fundamental question in biology, one dependent not only upon changes in species' environments (i.e. directional selection) but also on the amount of available phenotypic (co)variation. Our study investigated how these two interact over a period of approximately 100 generations in two co-distributed chipmunk species. We observed that some features of the available phenotypic (co)variation in cranial dimensions changed in response to directional selection, but in unpredictable ways.

Directional selection was a major component in skull evolution for both *T. alpinus* and *T. speciosus*, although the strength of selection, estimated as the selection gradient, for *T. alpinus* was stronger than for *T. speciosus*. The stronger selection gradient observed in *T. alpinus* populations supports the hypothesis that this species has been more sensitive to the environmental changes observed across its habitat range compared with *T. speciosus* [17,18,42]. As a species can become extinct when directional selection is too strong, one might think that *T. alpinus* is at a higher risk of extinction than is *T. speciosus*. Indeed, theoretical work has determined a threshold for the ratio of the amount of sustained environmental change, translated as the selection gradient, to the amount of variance available in a population above which the risk of extinction increases [43], and *T. alpinus* from both transects exhibited a ratio greater than this threshold [44]. However, as discussed below, available phenotypic (co)variation in *T. alpinus* has been redistributed between the historical and modern sampling periods to match the selection gradient, potentially enhancing the survival of this chipmunk.

The pattern of selection in each population is distinct. In *T*. *alpinus* from the Yosemite, most traits were selected to increase (electronic supplementary material, figure S2), which resulted in the observed size increase. On the other hand, *T*. *alpinus* from the Southern Sierras presented selection on most traits in the facial region to decrease, in congruence with the observed morphological changes in the population [22]. It is beyond the scope of our study to relate the selection gradients observed with possible selection agents. However, it is interesting to note that in all populations, one trait has been strongly selected to increase (IS–NSL; electronic supplementary material, figure S2). This trait (figure 2) represents the opening of the nasal cavity where the turbinate bones are located, which are responsible for water and heat regulation in mammals [45]. A possible explanation for this convergent change in all populations could involve selective pressure to increase the capacity to better regulate water and heat loss as climate has warmed. In Yosemite, for example, there has been as much as 3.5°C increase in average temperature, especially in the high country occupied by *T. alpinus*, over the past century [46].

The difference in selection strength between populations of *T. alpinus* and *T. speciosus* allowed us to further narrow our predictions about how the available phenotypic (co)variation is expected to change in a directional selection scenario under a model that allows for the evolution of the genotype–phenotype map [12,13,47]. As *T. alpinus* faced a stronger directional selection regime than did *T. speciosus*, we hypothesized that any changes in the available phenotypic (co)variation would be more pronounced in *T. alpinus*.

Our first prediction proposes that the amount of phenotypic variation would increase in the direction of the selective regime but not necessarily in other directions. Indeed, *T. alpinus* showed both increased variance in the direction of selection but not in other directions and an increase in the total amount of phenotypic variation for both populations examined. Furthermore, this increase was unique to each population because the direction of the net selection gradient is not the same on both *T. alpinus* populations (electronic supplementary material, figure S2). This was not the case for both populations of *T. speciosus*. In principle, increased (co)variation in the direction of selection is compatible with a hypothetical increase in the frequency of rare alleles, in which case it could be explained solely by additive genetic variance [10]. However, an increase in variance caused by rare alleles is thought to be transient and likely to most impact traits determined by a small number of alleles [48], which is unlikely to be the case for skull traits, and thus the impact of rare alleles is probably limited. Alternatively, a model accounting for epistatic interactions among traits could lead to this increased variation in the direction of selection, as indicated in figure 1. Under this model, selection acting in a given direction will favour alleles influencing the degree of correlation favoured by selection, which in turn will lead to changes in the amount of variation in this direction.

Our second prediction was that co-selected traits would increase their correlation. Indeed, Yosemite *T. alpinus* did show an increase in overall phenotypic integration among traits, confirming this prediction. This pattern is also in accordance with an epistatic model, where coordinate selection across multiple traits will lead to tighter correlations among them [6]. On the other hand, *T. alpinus* from the Southern Sierras did not exhibit an increase in overall correlation among traits. One possible explanation for these different spatial responses has to do with the direction in which selection acted, because selection in the Yosemite population was already in a dimension of relatively high variance, meaning that most traits were involved in the response to selection, while in the Southern Sierra population selection was not along a high variance dimension, as can be appreciated by comparing the evolvability boxplot distributions (figure 3). Therefore, we would expect fewer traits in the Southern Sierra transect to be co-selected, leading to the stability that we observed in the overall phenotypic integration among traits. Furthermore, Yosemite *T. speciosus* showed the opposite response, decreasing its overall degree of correlation, a pattern that would be expected in a drift scenario [11].

Lastly, we predicted a re-alignment of (co)variation patterns with the selection gradient. Even though this was not the case for three of the populations analysed (figure 5), the Southern Sierra sample of *T. alpinus* did exhibit increased correlation between PC1 and the selection gradient. Once again, this outcome might be linked to the direction in which selection was operating, because in the latter population, selection was not in a dimension of relatively high variation. On the other hand, for *T. alpinus* from Yosemite, the selection gradient was in a direction of high (co)variation, and thus a reorientation of the patterns of (co)variation would not have been necessary as sufficient variation was available for selection to act upon. These results are even more interesting if we link them to the evolvability results, as both *T. alpinus* populations increased the amount of variation along the selection gradient, but only one also re-aligned the matrix structure. Macroevolutionary studies in mammals have shown that the overall phenotypic correlation among traits and the amount of variance are very labile between groups, whereas the orientation of (co)variation is more conserved on a macroevolutionary time scale [24,25,49–52]. As we showed that a re-orientation of phenotypic (co)variation can be easily achieved under a model of sustained directional selection, a possible explanation for the widespread stability of the orientation of (co)variation on a macroevolutionary time scale could be a concordance between the adaptive landscape and patterns of (co)variation. On the other hand, this stasis pattern could be due to a high level of fluctuating selection, leading to the cancellation of directional selection effects on covariance matrices in the long run. An interesting next step would be to investigate additional cases where selection has not acted along an axis of major variance to see if the pattern reported here is robust.

An interesting aspect raised by our analysis is the striking contrast observed between species. While *T. alpinus* exhibited changes in the three aspects of its phenotypic (co)variation analysed, *T. speciosus* has remained fairly stable, with change observed only in the overall phenotypic correlation among traits, and in a direction opposite to the one we predicted. It is possible that the observed differences between these two species are related to the discrepant selection strengths they were subject to. Although we did observe a pattern consistent with directional selection for all populations, the selection pressure experienced by *T. speciosus* might have been weaker than that necessary to produce changes in phenotypic structure.

There are some caveats to our study that should be acknowledged. First, we worked with phenotypic instead of genetic (co)variances because of the difficulties in estimating the latter [51,53]. Although evolutionary quantitative genetics theory is based on genetic (co)variation, we assumed that our phenotypic estimates were appropriate substitutes for their genetic counterparts based on a substantial body of evidence showing that they are structurally similar, at least for morphological traits and particularly for mammals [23–25,54–57]. Furthermore, comparisons among the phenotypic (co)variation in the historical and modern samples of both chipmunk species and the genetic (co)variation of a third distant related rodent species (*Calomys expulsus*) showed that they were structurally similar, a result that supports our assumptions (for details on the reasoning behind this analysis, see [24]).

Second, we were able to estimate the net selection gradient only between the endpoints of the approximately 100 generations between the historic and modern periods. Although not ideal [58], this is the best approximation we have for the level of directional selection operating on both species between the two sampling periods.

Our study has several strengths. First, we examined well-sampled natural populations separated by multiple generations. Second, the large measured effective population size of our samples allowed us to overcome some of the caveats expected when working in experimental settings, which are frequently hampered by small effective sizes. Indeed, small effective sizes will affect any study designed to analyse the effects of directional selection because of the likelihood of substantial genetic drift. Moreover, the populations we examined were subject to differing selective regimes, in both direction and strength, over time. This allowed us to narrow the predictions and match our expectations to the different populations.

Our study examined a largely neglected aspect of evolutionary dynamics: the interaction between selective regimes and available phenotypic (co)variation. We suggest that the available multidimensional phenotypic (co)variation of a species can evolve quickly in natural populations under relatively strong directional selection, a hypothesis supported by both theoretical and simulation studies [12,13,59]. As species under strong directional selection tend to be more prone to extinction, our study, coupled with previous theoretical, computational and experimental knowledge, highlights one mechanism by which species may enhance their survival in the face of environmental change, namely a rapid reorganization of their available phenotypic (co)variation. This is especially relevant in the light of the ever-growing environmental threats to species survival.

## Data accessibility

Observed and estimated **P**-matrices and mean phenotypes for all populations are available on Dryad: http://dx.doi.org/10.5061/dryad.f8q6b [60].

## Authors' contributions

A.P.A.A., J.L.P. and G.M. designed the study. A.P.A.A. and G.M. delineated the analyses. A.P.A.A. measured the specimens and analysed the data. A.P.A.A., A.H., G.M. and J.L.P. wrote the manuscript. All authors read, commented on and approved the final version of the manuscript.

## Competing interests

The authors declare no competing financial interests.

## Funding

Fieldwork was financed by grants from the Yosemite Conservancy, The National Park Service Inventory and Monitoring Program, the National Geographic Society (Grant 8190-07) and the National Science Foundation (DEB 0640859). This work was supported by grants from the Fundação de Amparo à Pesquisa de São Paulo (FAPESP) to A.P.A.A. (2010/52369-0 and 2012/00852-4), A.H. (2012/24937-9) and G.M. (2011/14295-7).

## Acknowledgements

We thank Joseph Grinnell for his foresight in assembling the specimens and accompanying detailed data for the historical samples we had available, and members of the Grinnell Resurvey Project for the modern samples and data. We thank C. Patton, L. Chow, P. Moore and G. Lefèvre for fieldwork assistance; Miriam Zelditch, A. Kauffmann-Zeh, associate editor Rita Mehta and two anonymous reviewers for helpful comments on previous versions of the manuscript; G. Garcia for kindly providing the G-matrix; and D. Houle for helpful discussions.

## Footnotes

Electronic supplementary material is available online at http://dx.doi.org/10.6084/m9.figshare.c.3575741.

- Received July 19, 2016.
- Accepted October 28, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.