How does recent climate warming and climate variability alter fitness, phenotypic selection and evolution in natural populations? We combine biophysical, demographic and evolutionary models with recent climate data to address this question for the subalpine and alpine butterfly, Colias meadii, in the southern Rocky Mountains. We focus on predicting patterns of selection and evolution for a key thermoregulatory trait, melanin (solar absorptivity) on the posterior ventral hindwings, which affects patterns of body temperature, flight activity, adult and egg survival, and reproductive success in Colias. Both mean annual summer temperatures and thermal variability within summers have increased during the past 60 years at subalpine and alpine sites. At the subalpine site, predicted directional selection on wing absorptivity has shifted from generally positive (favouring increased wing melanin) to generally negative during the past 60 years, but there is substantial variation among years in the predicted magnitude and direction of selection and the optimal absorptivity. The predicted magnitude of directional selection at the alpine site declined during the past 60 years and varies substantially among years, but selection has generally been positive at this site. Predicted evolutionary responses to mean climate warming at the subalpine site since 1980 is small, because of the variability in selection and asymmetry of the fitness function. At both sites, the predicted effects of adaptive evolution on mean population fitness are much smaller than the fluctuations in mean fitness due to climate variability among years. Our analyses suggest that variation in climate within and among years may strongly limit evolutionary responses of ectotherms to mean climate warming in these habitats.
Numerous studies have documented biological responses to recent changes in climate, including changes in seasonal timing, migration, phenotypic plasticity, geographic distribution and range edges, population abundances and habitat use . For organisms with limited potential for migration, evolution may also enable important responses that allow populations to adapt to local and changing climates. Evolutionary responses to recent climate change have been documented in several systems, for phenotypic traits that include seasonal timing [2–4], coloration  and thermal sensitivity [6,7]. These studies demonstrate that organisms can evolve in response to climate warming over the timescale of one to a few decades. Conversely, other historical studies have failed to detect evolutionary responses to climate change [1,8]. One limitation is that such historical studies do not provide a quantitative, mechanistic understanding of how specific aspects of climate variation and change affect the dynamics of phenotypic selection and microevolution of traits over time. In addition, the fitness consequences of such evolution are poorly understood: how much does adaptation to environmental change improve mean population fitness ?
A theoretical framework for describing evolutionary responses to sustained, directional environmental change (including climate change) is now well established [10,11]. These models consider stabilizing selection towards some optimal trait value, where the optimal trait value changes directionally with time. In this scenario, the mean trait value lags behind the optimal value, generating selection and (potentially) evolution of the population; if the lag becomes too great, the population will go extinct. The degree of evolutionary adaptation and the probability of extinction depend on the rate of environment change, genetic and phenotypic variation, and environmental and demographic stochasticity . These general models have provided important qualitative insights into the factors that influence adaptation and extinction in response to climate warming. However, predicting these responses for specific empirical systems requires a quantitative understanding of how climate determines optimal trait values, fitness functions and selection. The few applications of these models to empirical systems demonstrate that trait evolution can maintain fitness and distributions during climate change .
Climate change is expected to increase environmental variation in addition to shifting mean environmental conditions . This variation, both within and across generations, drives fluctuations in selection [14,15], which can slow evolution in response to directional environmental change. This effect is complicated by fitness integrating environmental responses at multiple timescales. Fitness of our focal thermoregulatory trait—wing melanin—integrates responses to chronic environmental conditions (via time available for flight) and acute thermal stress events (via overheating and declines in egg viability). Geographic differences in wing melanin are important for thermal adaptation to climate in a number of butterfly species [16–19].
In this study, we explore selection and evolution of wing melanin in the alpine and subalpine butterfly, Colias meadii, in the southern Rocky Mountains (USA). Many species in these habitats are adapted to local climate conditions and have limited potential for large-scale migration and gene flow. We integrate microclimate, biophysical and demographic models to map phenotype onto fitness in different and variable climatic conditions. We combine these with a simple quantitative genetic model to predict selection and evolutionary responses to climate change. Our analyses have three main goals. First, we predict fitness functions, selection and evolutionary responses in wing melanin in relation to climate changes in this region during the past 50–60 years. Second, we compare the temporal patterns of selection and evolution at subalpine and alpine sites. Third, we evaluate how population mean fitness at these sites changes with time, and whether evolutionary adaptation can increase mean fitness or reduce variance in mean fitness over time. Our results suggest that climate changes can generate substantial directional selection on this trait, but that stochastic variability in climate strongly limits the evolutionary responses to this selection.
2. Material and methods
(a) Study system and thermal biology
We consider C. meadii, a species native to the mountains of western North America. In the southern Rocky Mountains, C. meadii occurs in alpine and subalpine habitats at elevations between 2900 and 4300 m. It has an obligate larval diapause and a single adult flight season in July–August each year. Here, we focus on subalpine and alpine populations in the lower and middle parts of its elevational range in Colorado. We chose this species in part because of the availability of phenotypic data on key thermoregulatory traits collected in 1980 . We emphasize that our analyses are intended to be representative of the responses of this species in subalpine and alpine regions in Colorado, rather than to predict responses at specific locales.
Colias are strong flyers, and active flight is essential for courtship and mating, nectaring, oviposition and other activities [21–23]. Flight activity is restricted to body temperatures between 27 and 40°C and flight performance is greatest at 33–38°C [16,20] (H. MacLean 2014, unpublished data). The thermal sensitivity of flight is similar for different populations and species of Colias. Adults behaviourally thermoregulate to achieve the body temperatures needed for flight, and do not use endogenous heat production to elevate body temperatures . Colias use a lateral basking posture with the wings closed and the ventral hindwing (VHW) surfaces oriented perpendicular to the sun to increase body temperatures. At body temperatures above 40°C, butterflies stop flying and use a heat-avoidance posture with the body and wings oriented parallel to the sun; they also crawl down into the vegetation to avoid direct solar radiation.
Colias populations and species are adapted to local climates through two key morphological, thermoregulatory traits [16,24,25]: the degree of melanism on the VHWs (determining solar absorptivity); and the length of setae on the thorax (fur thickness). Wing absorptivity is determined by the relative proportions of melanic (black) and pteridine (yellow or orange) scales . Colias populations and species from higher elevations and latitudes have greater wing solar absorptivity and thoracic fur thickness, enabling them to achieve the elevated body temperatures needed for flight even under cool environmental conditions [20,27]. Extreme high temperature can also affect survival and fecundity. For example, a daily heat shock at 45°C significantly reduced adult lifespan and caused a fourfold reduction in egg production . Importantly, Colias adults may be exposed to short intervals of deleteriously high body temperatures (more than 40°C) even at high elevation where temperature variation can be substantial .
(b) Climate and microclimate
We used daily maximum and minimum temperature data from two weather stations, to represent climate conditions at the lower (subalpine) and middle (alpine) parts of the elevational range of C. meadii in Colorado. Data for the subalpine site (site C1: 3048 m, 40.03 N, 105.55 W) derive from the Niwot Ridge LTER (http://niwot.colorado.edu) that were cleaned and filled . Data for the alpine site (Climax: 3514 m, 39.38 N, 106.20 W) derive from the National Weather Service Cooperative (COOP) Program. We estimated air temperatures (Ta) at 10 min intervals using a diurnal temperature variation function incorporating sine and exponential components . For simplicity, we refer to C1 as the subalpine site, and Climax as the alpine site.
Global horizontal solar radiation was calculated as a function of elevation, latitude and longitude by discounting global extraterrestrial radiation . Radiation was then partitioned into direct and diffuse components as a function of the atmospheric transmissivity τ (ratio of global horizontal solar radiation at surface and calculated global extraterrestrial (top of atmosphere) horizontal solar radiation). Tau distributions were estimated hourly using several years of data from the NREL Solar Radiation Research Laboratory Baseline Measurement System in Golden, CO, USA (1829 m, 39.74 N, 105.18 W http://www.nrel.gov/midc/srrl_bms/). We used kernel density estimation to simulate a τ value for each time interval. Solar radiation was partitioned using an empirical relationship  as modified for high-altitude sites in Colorado .
We implemented a microclimate model [34–36] using finite-difference methods to solve heat balance equations describing soil temperatures at the surface and specified depths. We scaled microclimate variables to plant heights of 0.2 and 0.5 m by estimating temperature and windspeed profiles  using data collected at heights spanning 0.05–1.5 m in July 2012 at the subalpine site (see electronic supplementary material). Based on weather station data from July 2011 at this site, the mean windspeed at 0.5 m height was 0.4 m s−1. We used this value in our models for both sites. Because alpine sites can have higher windspeeds than subalpine sites (depending on local topography), our model predictions may underestimate optimal wing absorptivities at the alpine site (see Results).
(c) Biophysical and demographic models: connecting climate, traits and fitness
We use a steady-state heat flux model for Colias that was developed and field validated by Kingsolver [20,21] to predict thoracic body temperature (operative environmental temperature, Te) based on thermoregulatory traits (body size, basal ventral hind wing solar absorptivity and thoracic fur thickness), behavioural posture (basking and heat-avoidance) and environmental conditions . The range of possible values of wing absorptivity α in Colias is from 0.4 (all yellow pteridine scales) to 0.7 (all black melanic scales). For both subalpine (C1) and alpine (Climax) sites, we assume a fur thickness of 1.46 mm and thorax diameter of 3.6 mm, based on measurements for C. meadii at several sites in Colorado . We assume that butterflies select the body temperature closest to their thermal optima (35°C) with available body temperatures bracketed by those in full sun and full shade (no direct radiation). The model successfully predicts patterns of Te, flight activity time and heat-avoidance in the field for C. meadii and other Colias species along an elevational gradient in Colorado [20,38]. Predictions of Te are updated every 10 min.
Buckley & Kingsolver  used this biophysical model together with a demographic model to connect climate and thermoregulatory traits to predictions of fitness (net reproductive rate, λ) (see electronic supplementary material for details). We calculate daily egg production per female as the product of available flight time and the rate of oviposition. We estimated the probability of flight as a plateauing function of operative environmental temperatures, Te: Pflight = exp(−0.5 × (abs(Te−33.5))/5)3.5. The function is based on flight data for Colias eriphyle in Montrose, CO . Lab experiments show that high body temperatures can reduce the viability of oocytes in female Colias . We incorporate egg viability as a function of body temperature, modelled as an exponential decay from 1 at 40°C to 0.75 at 50°C . We multiply daily egg production by the average of hourly viability estimates. We estimate λ by summing over days to either a duration of 5 days, reflecting the mean adult life span in the field [40,41], or reaching a maximum lifetime egg production of 700  as the product of survival to maturity, daily survival and egg production. We simulated 500 individuals with a uniform probability of initiating adulthood across a single annual generation with a single summer flight season .
(d) Selection and evolution models
Our microclimate, biophysical and demographic models allow us to predict the fitness λ of an individual Colias as a function of climate variables and thermoregulatory traits. We use a simple quantitative genetic model to predict selection and evolution of solar absorptivity α of the VHW in response to climate. We used estimates of the phenotypic mean of α for C. meadii (0.649) reported by Kingsolver [20,21] for population samples taken in 1980. Ellers & Boggs  used parent–offspring breeding experiments to estimate the narrow-sense heritability h2 of α for C. eriphyle, yielding h2 = 0.43 for males and 0.36 for females ; we use a value of 0.40 in our simulations. We assume that selection is sufficiently weak that the heritability and phenotypic and genetic variances do not change with time . We used more extensive recent sampling to estimate the standard deviation of α as 0.062 (H. MacLean 2014, unpublished results). Additional simulations (not shown) demonstrate that the precise values of h2 or phenotypic variance have little effect on our quantitative or qualitative results. We also assume that thermal conditions do not produce phenotypically plastic changes in α (see Discussion).
Our models predict the fitness function relating solar absorptivity α to fitness for a given site and year. Combined with estimates of the phenotypic mean and variance, we can estimate the (unstandardized) directional selection gradient β for α, and use the heritability h2 to predict the evolutionary response to selection (change in mean phenotype) in the next generation .
(a) Changes in climate
The monthly mean of daily maximum temperatures in July has increased significantly at both sites (0.6 and 0.4°C/decade, for the subalpine and alpine sites, respectively) during the past 55 years (figure 1a), in agreement with previous climatic analyses in this region . The variation in temperature among years at each site is striking. Indeed, the total increase in mean daily maximum temperatures over the past 50 years (3 and 2°C, for the subalpine and alpine sites, respectively) from 1955 to 2010 at each site is similar in magnitude to the standard deviation among years around that warming trend (1.57 and 1.28°C, respectively). The within-year variation (standard deviation) in July daily maximum temperature has also increased significantly during the past 60 years (figure 1b), again with substantial variation among years; and the within-year standard deviation is larger than the mean warming at each site and is particularly pronounced at the subalpine site. As a result, we must consider climate warming at these sites in the context of the substantial temporal variation in climate within and among years.
(b) Fitness functions
At the subalpine site (figure 2a), fitness functions have generally shifted towards the left between 1955 and 2010, towards lower optimal α and higher maximal fitnesses. However, there is substantial year-to-year variation in the position and height of the fitness functions (figure 2a). Fitness functions at the alpine site show a different pattern (figure 2b). In most years, fitness at the alpine site increased monotonically with absorptivity, yielding an optimal α of 0.70 (figure 2b). Conversely, in many years, values of α below 0.55 yielded fitness values below replacement rate (i.e. λ < 1.0). As a result, increased wing melanin has generally been favoured at this site during the past 60 years.
(c) Selection, evolution and mean fitness
We next use these fitness functions to predict temporal patterns of directional selection and evolutionary responses (see Material and methods). At the subalpine site (figure 3a,b), the directional selection gradient β has declined substantially. The selection gradient for α was consistently positive (selection for increased wing melanin) from 1955 to 1975, but was negative (selection for decreased wing melanin) during most years from 1987 to 2010 (figure 3b). We predicted little evolutionary change in mean α during the period from 1980 to 2010, largely because of the fluctuations in the direction and magnitude of directional selection (figure 3a,b).
The predicted directional selection gradient for α also declined at the alpine site (figure 3c,d). However, directional selection was always positive (selection for increased wing melanin), reflecting the fact that the optimal α was 0.7 for nearly all years during this period. The magnitude of selection varied strongly among years, but was generally greater at the alpine (figure 3c) than at the subalpine (figure 3a) site. Our 1980 empirical estimate of α below the predicted optima resulted in the prediction that α would increase rapidly until 2010, as a result of the strong directional selection.
(d) Fitness consequences of climate change and evolution
Both the subalpine and alpine sites had substantial variability in predicted mean fitness among years (figure 4). In general, overall mean fitness was 6–8% higher and variation (standard deviation) in mean fitness was 37–48% lower at the subalpine compared with the alpine site. Mean fitness increased significantly with time at the alpine site (F = 8.150, p = 0.006) (figure 4a), but not at the subalpine site (F = 2.154, p = 0.148) (figure 4b). To highlight the potential effects of evolutionary adaptation on mean fitness, we compare two scenarios that differ in the heritability of absorptivity, one without evolution (heritability h2 = 0) and one with maximal evolution (h2 = 1) (figure 4). The evolution scenario did yield increased overall mean fitness and reduced among-year variation in mean fitness at both sites. Relative to no-evolution, the evolution scenario increased the overall mean fitness by 1% (subalpine site) and 3% (alpine site). Evolution had greater effects on variation in mean fitness: evolution reduced the standard deviation in mean fitness by 23% (subalpine site) and 7% (alpine site). The predicted effects of evolution on the overall mean and variation in mean fitness are considerably smaller than the variation in fitness caused by variation in climate among years (figure 4).
Understanding how environmental variability alters selection and evolutionary responses to sustained, directional environment change is an important challenge for evolutionary ecology. Our analyses illustrate the benefits of linking variation in climate, phenotype and fitness in addressing this challenge for empirical systems. As expected, mean climate warming during the past half century at subalpine and alpine sites has altered predicted patterns of directional selection and reduced optimal absorptivities, but these mean changes are similar in magnitude to the stochastic variation in selection and optimal trait values among years. Climate warming also generates complex changes in the fitness functions beyond reductions in the optimal trait value (figure 2). For example, lower optimal absorptivity is associated with higher maximal fitness at the optimum (r = −0.36 at the subalpine site). This pattern, consistent with the ‘hotter is better’ hypothesis that populations adapted to higher temperatures may achieve greater maximal fitness [45,46], highlights that warming may differentially affect absolute and relative fitness. In addition, lower optimal absorptivity is associated with higher (less negative) curvature of the fitness function (R = −0.48 at the subalpine site). The curvature of the fitness function results from two opposing components: darker wings enable additional flight and reproductive output, but increase the risk of thermal stress that reduces egg viability and adult survival. Our analyses emphasize how these two components interact with environmental variability in ways that would not be anticipated using simple theoretical models of evolutionary responses to directional environmental change.
A key outcome of our study is that, despite substantial mean climate warming at these sites, predicted evolutionary responses of wing absorptivity to warming are modest. Multiple factors contribute to this result. First, stochastic temporal variation in climate and selection will reduce the evolutionary response and increase the lag between the mean and optimal absorptivity, as predicted by theoretical models . Note that thermal variability within years at these sites is also increasing with time (figure 1), further reducing evolutionary responses. Second, the negative associations of optimal value of absorptivity with maximal fitness and the curvature of the fitness function will reduce the strength of directional selection on the trait. Third, the fitness function is generally asymmetric, with a shallower slope above than below the optimal trait value (figure 2). As a result, directional selection near the fitness peak will be stronger on absorptivities below than above the optimum. All of the factors will tend to reduce the evolutionary response to mean climate warming during the past half century in this study system, and potentially increase the evolutionary lag between mean and optimal trait values. How these evolutionary responses and lags are likely to change in response to future climate warming remains unexplored.
Climate variation also had important and unexpected effects on predicted variation in mean population fitness. Temperature variation within and between years was generally larger (figure 1), but variation in mean fitness was generally smaller (figure 4), at the subalpine than the alpine site. This nonlinear relationship between temperature and fitness variation is not due to heat stress at higher temperatures, as reported in some previous analyses of climate change in ectotherms [47,48]. Instead, relatively small increases in temperature at cool, alpine sites can have large effects on increasing opportunities for flight and reproduction. Unlike at the subalpine site, recent warming has ameliorated the climate conditions for adult Colias in alpine habitats , as reflected in increases in predicted mean fitness over the past 60 years.
At the alpine site, we find consistent, but declining, positive directional selection leading to evolutionary increases in mean wing absorptivity during the past 30 years (1980–2010). Our prediction of high values for optimal absorptivity may stem from incompletely accounting for the influence of thermal variation on body temperatures. For Colias butterflies and many other large insects, changes in solar radiation on the timescale of minutes to hours can directly impact body temperature and performance [28,50]. Limited past data for cloudiness and radiation and the uncertainties of future scenarios are major limitations on our ability to predict the consequences of climate change for terrestrial ectotherms. Our models may underestimate the extent to which spikes in solar radiation result in thermal stress.
Our models for selection and evolution in adult butterflies do not consider the potential effects of climate change on other life stages. Climate change can clearly have important ecological effects on multiple life stages in insects [1,51,52]. Additionally, a recent study documented evolutionary changes in thermal sensitivity for larval feeding in other Colias species in response to recent climate changes . Integrating the potential effects of climate change across multiple life stages remains an important challenge .
Our current models do not consider the potential effects of phenotypic plasticity on trait values, selection and evolution. For some Colias species with multiple generations per year (multivoltine), day length or temperature during late larval and pupal development can influence adult wing melanin, resulting in different seasonal phenotypes [25,54,55]. Plasticity of wing melanin has not been studied in C. meadii and other univoltine species of Colias that have an obligate winter diapause. Theoretical models increasingly recognize the potential for phenotypic plasticity in response to environmental variability to modify evolutionary responses . Adaptive plasticity may initially ameliorate the negative fitness consequences of environmental change, but may also reduce rates of evolutionary responses and increase the evolutionary lag of the population. Although we do not anticipate that considering adaptive plasticity would alter our results qualitatively for univoltine Colias in subalpine and alpine habitats, incorporating empirical data about plasticity into models that link phenotypes to fitness in this or other systems will help resolve the role of plasticity in responding to environmental change.
The potential for evolutionary adaptation to climate change to improve the mean fitness of populations and thus reduce the risks of extinction remains an important issue for evolutionary ecologists . Our results for Colias butterflies suggest that evolution does have the potential to reduce variation in mean fitness, but they predict that evolution has had little impact on increasing overall mean fitness during the past half century in these systems. In general, among-year variation in mean fitness due to climate variability was substantially greater than the effects of evolution in subalpine and alpine Colias. The role of stochastic, fine-grained variation in climate in producing plastic responses to climate change has been widely documented and discussed, but the impact of such variation in reducing evolutionary responses to mean climate warming has received less attention [8,56,57]. Numerous analyses suggesting that potential rates of evolution are too slow relative to rates of climate change largely omit consideration of variability in climate. Such considerations are likely to further erode the potential for evolution to ameliorate responses to climate warming in many systems. The potential fitness consequences of evolution in response to future climate warming must be considered in the light of variability in climate.
Research supported in part by NSF grant no. DEB-1120062 to the authors.
We thank Jessica Higgins, Ray Huey, Heidi MacLean and Rory Telemeco for useful discussion and comments on the manuscript, and the Rocky Mountain Biological Labs for research facilities during part of the research.
- Received October 8, 2014.
- Accepted December 23, 2014.
- © 2015 The Author(s) Published by the Royal Society. All rights reserved.