Historical changes in thermoregulatory traits of alpine butterflies reveal complex ecological and evolutionary responses to recent climate change

Trait evolution and plasticity are expected to interactively influence responses to climate change, but rapid changes in and increased variability of temperature may limit evolutionary responses. We use historical specimens to document changes in the size and thermoregulatory traits of a montane butterfly, Colias meadii, in Colorado, USA over the past 60 years (1953–2012). We quantify forewing wing length, ventral wing melanin that increases solar absorption, and the length of thorax setae that reduces convective heat loss. The mean of all three traits has increased during this time period despite climate warming. Phenological shifts may have extended the active season earlier at low elevations and later at high elevations, increasing exposure to cool temperatures and selecting for increases in thermoregulatory traits. Fitness increases at higher elevations due to warming could also increase thermoregulatory traits. Warmer temperatures during pupal development and later flight dates in the season are associated with decreased wing melanin, indicating a role of phenotypic plasticity in historical trait changes. Phenotypic shifts result from a complex interplay of ecological and evolutionary responses to climate change. Environmental variability within and across seasons can limit the evolutionary responses of populations to increasing mean temperatures during climate change.


Background
For many organisms, the rate of climate change may be sufficiently fast to preclude evolutionary adaptation [1,2]. However, there is some evidence for evolutionary changes in body size [3], coloration relative to background matching [4][5][6], and phenology; specifically the timing of diapause [7], of nesting [8] and of flowering [9]. Evidence for evolutionary or historical changes in traits relevant to thermal physiology is limited (but see [10][11][12][13]) Natural history collections offer a largely untapped resource to examine morphological responses to climate change [14]. In particular, insects are abundant in these collections and the thermal significance of morphological traits has been extensively documented for some species.
A biophysical model for Colias butterflies enables estimating how morphological thermoregulatory traits and environmental conditions interact to determine body temperatures [15][16][17][18]. Specifically, the body temperatures of butterflies increase proportionally to increases in air temperature [17]. Two key thermoregulatory traits influence body temperatures: 1) Increased melanism on the posterior ventral hindwings allows the butterflies to absorb more solar radiation via closed wing basking, and 2) longer setae on the thorax reduce convective heat loss [19,20]. Wing melanin is heritable and increases with elevation for Colias species [21], suggesting its importance in determining fitness. Moreover, the biophysical model enables linking phenotypes to function because Colias flight is restricted to a narrow range of body temperatures (28°-40°C) [17,20]. Thus their flight time is often limited by climate, especially at higher elevations and latitudes [22][23][24]. Flight time is a determinant of reproductive success because flight is required for activities including mating, feeding, and oviposition [16,19]. Conversely, butterflies may overheat, which not only limits flight time, but can also directly reduce the number of viable eggs produced by each female [25]. Past Colias research provides a mechanistic understanding of how environmental conditions interact with thermoregulatory traits, which are well preserved in collections, to determine fitness.
Application of the biophysical model suggests that recent warming in Colorado has increased both flight time and overheating risk for a montane Colias species, C. meadii [16]. Translating estimates of flight time and overheating risk into fitness based on egg-laying suggests that recent warming has driven declines in fitness at lower elevations and fitness increases at higher elevations for C. meadii with constant thermoregulatory traits [16]. An analysis of selection gradients complicated our initial expectation that climate warming would select for decreased melanism. The analysis estimated that directional selection for decreased wing melanism has become more consistent and stronger in the lower elevation subalpine portion (3.04 km site) of C. meadii's distribution over the last 60 years in Colorado [26]. The analysis predicted consistent selection for increased wing melanism at the high elevation extent (3.51 km site) where flight time is more limited, but the prediction could be an artifact of insufficiently capturing temperature extremes that could decrease survival of individuals with high wing melanism. The model also indicates that seasonal and annual variation in climate generates substantial variation in the strength and direction of predicted selection, which would dampen the potential for phenotypic shifts in response to recent climate change [26]. To address these predictions, we measured the ventral hindwing melanism, thorax setae length and wing length (proxy for body size) in historical specimens. Available specimens were restricted to a fairly narrow elevation range (25 to 75% quantile: 3.36 km to 3.64 km) spanning the high elevation site in the model. Thus, we expect that the direction of selection may vary with elevation across our sampled populations and that environmental variability may act to dampen directional selection.
Phenotypic plasticity is also likely to play an important role in responding to climate change, particularly in environmentally variable regions [2,27]. Plasticity may either facilitate evolution through enabling persistence or hinder evolution through buffering selection [28]. Montane Colias exhibit plasticity where increased temperature during pupal development decreases melanism [10,29]. Thus, we expect directional selection and phenotypic plasticity to interact for montane Colias.
Climate warming is also expected to shift insect body size via selection and plasticity. Warmer temperatures accelerate development, which can shorten the duration of development and reduce adult size (temperature-size rule) [30]. Warmer environments can also select for decreased size to decrease body temperatures [17]. Consequently, predictions and observations of insect responses to climate change generally indicate a decrease in body size [3,31]. However, spring warming can extend the time available for development and increase size [32].
Here we examine the roles of evolution and plasticity in responding to climate change by using museum specimen collections to document shifts over the past 60 years (1953-2013) in morphological traits in C. meadii butterflies in Colorado. We evaluate several factors expected to determine trait values via both phenotypic plasticity and selection. We examine developmental temperatures, which determine phenotypic plasticity, and temperatures during the previous flight period, which determine the reproductive success of parents. We expect mean trait values to decrease with the calendar date of collection because mean temperatures increase and daily variability decreases through the season (analysis of data from http://www.wrcc.dri.edu/climatedata/climsum accessed 4 May 2015). Given the seasonal progression of temperatures, we also consider the interactions between developmental (or pupal) temperatures and collection date. Finally, we include collection year as an indicator of evolutionary shifts in phenotypes. However, we are unable to control for other factors that vary over long time scales such as trends in temperature mean and variability omitted from our temperature measures; trends in habitat, host plant availability, wind speed, or cloud cover; and collection biases.

Trait measurements
Colias meadii occurs in subalpine and alpine meadows above 2.50 km elevation throughout the Rocky Mountains, from northern New Mexico, USA to southern Alberta, Canada [20]. They have a single generation each year (univoltine) and diapause over the winter as 3 rd instar larvae; the adult flight season lasts 4 to 6 weeks in July and August (between Julian days 180-240). The subspecies C. meadii meadii is found largely in the southern part of the range at elevations between 3 and 4 km. It is genetically distinct from more northern populations of the species in Wyoming Basin [33]. In this study we focus on C. m. meadii in the southern part of the range, where summer temperatures have increased significantly over the last 60 years [16]. In this region, C. m. meadii individuals pupate in early to mid-June (between Julian days 150-170), and adult butterflies emerge in early July.
We processed specimens from museums in the United States with the largest holdings of C. m. meadii from Colorado (the Yale Peabody and the McGuire Center for Lepidoptera and Biodiversity); museums located in Colorado (the C.P. Gillette Museum of Arthropod Diversity at Colorado State University and the University of Colorado Museum); and museums with holdings from time periods absent from other museums (the Smithsonian Museum of Natural History, the Milwaukee Public Museum, and the California Academy of Science). Because museum holdings of C. m. meadii are sparse for the past 25 years, we also solicited samples from private collectors via the Yale Lepidoptera list serve and measured samples from one private collector.
Because the collections contained few females after 1990, we restricted our analyses to males (N = 379 from 20 counties). While females tend to have darker wings than males, elevation clines in melanism are similar between the sexes [21]. Model predictions for temporal trends in melanism are based on the temperature dependence of flight time and egg viability, which were translated into female reproduction [16,28]. The flight time constraints should influence male fitness similarly, so predicted temporal trends also apply to males. The similar elevation clines between sexes suggest that they respond to climatic gradients similarly and our observations should be indicative of both sexes.
All specimens for the C.m.meadii range were included except two specimens with a Julian collection date of 170 (June 19, collected at Como, Colorado 39.31 N and 105.89 W in 1984, McGuire Center collection) and 171 (June 20, collected at Los Pinos Pass, Colorado 38.1 N and 106.97 W in 1967, Yale Peabody collection), because they strongly influenced regressions and the date was more than two standard deviations below the mean (mean ± sd = 208.4 ± 13.65). Including these specimens increased the slope of our regressions and significance of our results.
Pinned specimens were measured by removing labels with forceps, transcribing the locality data, and then placing the head of the pin in a lump of modeling clay to expose the ventral hindwing (the wing region accounting for the majority of absorption for this ventral basking species). Elevation was recorded from the locality tag in the collections if available; otherwise we used a digital elevation model to estimate elevation based on latitude and longitude. We used Google Earth to georeference locality descriptions when necessary.
A digital micrometer was used to measure the forewing of the specimen from the thoracic insertion to the apex of the wing. The specimen was then photographed through a 100-mm macro lens in RAW format with a Canon Rebel XSi mounted on a copy-stand. Each image included a black and white standard. Because the height of the animal on the pin was variable, we used auto focusing to allow for the clearest image. We measured setal length on the ventral thorax with an ocular micrometer on a Wild M5 microscope as the longest setae between the first and second leg. All measurements were taken by MacLean and specimens were prepared for measurements in groups of five to obscure each specimen's metadata during measurements. Specimens were distributed randomly throughout collection boxes and thus were measured randomly within each collection.
We photographically assessed the degree of wing melanism on the posterior ventral hind-wing. First, we selected a triangulated region between the eyespot, hind wing insertion, and the wing margin [21,34]. Using a MatLab program (T. Hedrick, unpublished), we then converted the RAW image to black and white to account for potential fading. Black and white is appropriate because wing scales are pigmented by either melanin (dark) or pterin (light). We digitized the region of interest and the black and white standards for each sample, and then used the standards to calculate a standardized grey-level value between 0 (white) and 1 (black). The grey-level value should be robust to butterfly size because it represents the proportion of melanic wing scales. In order to verify our measure of grey-level as a proxy for absorptivity, samples collected in the field 2012 and 2013 were photographed, and the absorption spectrum was measured from 350 to 1050 nm in a spectroreflectometer with an optical integrating sphere (FieldSpecPro FEFR 7501, ASD Inc) for the same wing region. Previous studies used absorption at 650 nm as a measure of melanism and solar absorptivity for Colias wings (Watt 1968). Grey level and absorption at 650 nm were highly correlated (n = 60, R 2 = 0.78), confirming that grey level is an appropriate measure of solar absorptivity.

Climate data
To evaluate temporal changes in climate, we selected a representative NOAA Cooperative weather station (data retrieved from the Western Regional Climate Center, http://www.wrcc.dri.edu on 4 May 2015) at 3.5 km in Climax, CO (39.37 N, 106.18 W). This weather station was selected for its relevant elevation, its central location among sampling sites (median distance 61 km, 10 th quanitile = 23.3 km and 90 th quanitile = 96.0 km), and the completeness of its climate records throughout the time period of interest. We chose to use temperature data from one representative weather station because long-term data are lacking for most stations above 3 km in Colorado and, while interpolated data sets offered more complete coverage, they can perform poorly at capturing temporal variation for high elevation sites in the Rocky Mountains [35].

Data analysis
We estimated temperatures expected to influence each focal trait based on lab observations for the closely related Colias eriphyle. Melanism (and we assume setae length) is influenced by pupal temperatures (F 1,54 = 6.4, p < 0.01 [36]) and pupation lasts for 14 days at 20°C. We thus estimated temperatures during pupal development using data from 26 days to 6 days prior to collection (with the slightly extended duration intended to account for uncertainties). Forewing length is influenced by temperature during both late larval development (determining pupal mass: F 1,108 = 676.5, p < 0.0001) [36]) and pupation (determining wing length: F 1,101 = 6.5, p = 0.01 [36])). We thus estimate developmental temperatures using data from 46 days to 6 days prior to collection. We consider temperatures experienced by adults during the previous flight season, which determines reproductive success via flight time and overheating. We estimated temperatures during the previous flight season as the mean across Julian days spanning 201 to 216. The flight season was determined by taking the 25 th to 75 th quantiles of collection date. We used a constant flight season across years because we had little basis for predicting when in the flight season a parent flew. We examined sensitivity to our selection of time windows for developmental and pupal temperatures and found the model results to be robust. We additionally considered substituting flight season temperatures from the previous year with temperatures during the current year. Temperatures from the previous year were stronger predictors, but including current year temperatures did not qualitatively alter our model. We focus on temperature, but note that shifts in cloudiness and radiation could also influence plasticity and evolution. Increases in cloudiness on the order of 1.4% (of the sky) per decade have been estimated for the United States between 1976 and 2004 [37]. However, shifts are regionally variable. Analysis of solar radiation data proximate to our specimens suggests that solar radiation has increased over our observation period [38], consistent with other observations for the Rockies [39]. However, all these estimates have very low confidence due to a shift in techniques for cloud observations (which also influences solar radiation data derived from modeling) [37]. The radiation trends are sufficiently uncertain and of a magnitude that justifies our focus on temperature. Additionally, our focus does not include wind speed. While wind plays a role in the heat budget for these butterflies, there is little evidence (beyond a weak change in mean at high elevation) of a temporial trend at these sites [40].
We account for uneven collection intensities by restricting our analysis to 15 individuals per site per year. We then averaged the model output to perform model selection. We randomly sampled specimens each of the 50 times we repeated the analysis. Statistical model results were similar across iterations of the analysis, so we report results averaged across the 50 iterations. We performed model selection and averaging using the R package MuMIn [41] (commands: dredge and model.avg). For each iteration of the analysis, we selected the top 20 sub-models according to AICc (sample size corrected Akaike Information Criterion). All predictor terms were included among the best (delta AICc < 2) sub-models and the full model received strong support (AICc < 0.3). As a result, we focus our results on the full model (trait~developmental T (or Pupal T) + previous flight season T + Date of Collection + Year + developmental T (or Pupal T) x Date, where T is temperature in°C), but additionally report the percent of the top sub-models that included each predictor term. We present the full results of the model-averaging in Additional file 1: Table S1. We normalized the predictor variables using the scale function in R. We assessed collinearity using the variance inflation factor (vif function from the R car package). We did not detect problematic collinearity among our predictor variables (vif <2 for all main effects). Model selection excluded elevation as a poor predictor of trait values, likely because collections were concentrated over a fairly narrow elevation range (25 to 75% quantile: 3.36 km to 3.64 km).
To account for trait similarity due to geographic proximity, we used maximum-likelihood spatial autoregressive (SAR) models ( [42], R package spdep [43]). We selected an SAR error model over a lag or mixed model following assessment using a Moran's I test on the residuals of linear models as well as a Lagrange multiplier specification test. We used Moran's I tests, estimates of model performance, and spatial correlograms to select a 40 km threshold distance when developing spatial neighborhoods based on latitude and longitude. Neighbors were weighted using row standardization. We used a Moran's I test on the model residuals to confirm that the SAR models fully accounted for spatial autocorrelation in the data [42]. We report log-likelihood tests comparing the SAR to a null model consisting only of an intercept along with Nagelkerke pseudo R 2 values. We compared the SAR model output to linear models and found no qualitative differences in the direction of the effects.
We used an ANOVA to evaluate phenological trends over time. We standardized both year and elevation to the mean. Model selection identified elevation as an important predictor of phenology, despite its being omitted from the best models for trait values.
Influences of temperature and collection date on thermoregulatory traits suggest plastic responses to environmental temperatures. Non-significant tendencies for melanism to decrease with increases in pupal and flight season temperatures (Fig. 2) have high relative importance across models (relative importance values of 0.66 and 0.99, respectively: see Additional file 1: Table S1). Melanism decreases significantly with collection date, such that lower melanism is found later in the season (relative importance of 0.92). An interaction between pupal temperatures and collection dates also receives moderate relative importance across models (relative importance of 0.38). The interaction represents a correction accounting for the tendency for butterflies collected later in the season to have experienced warmer temperatures during pupation. Developmental temperatures and temperatures during the previous flight season received support as predictors of wing length (relative importance 0.89 and 0.34, respectively) with a tendency for warmer temperatures to increase size (Fig. 1). The interaction between developmental temperatures and collection date received some support as a predictor of forewing length (relative importance of 0.27). Setae length is increasing in response to increased temperatures in the previous flight season ( Table 2). Neither pupal temperatures nor their interaction with collection date are significant predictors of setae length, but they receive support for model inclusion (relative importance of 0.85 and 0.24, respectively, Fig. 3).

Discussion
Phenotypic plasticity is increasingly viewed as a primary response to environmental change [2,44,27]. Plasticity is increasingly being incorporated in predictive models of phenotypic change and warming tolerance [45] and being explicitly tested for in empirical studies [2]. The interaction between phenotypic plasticity and adaptive evolution in responding to environmental change remains unclear: does phenotypic plasticity tend to facilitate adaptive evolution by enabling persistence or does it tend to slow adaptive evolution by buffering selection [2,45]?  Fig. 1 Partial residuals (residuals of regressing the response variable on the independent variables, but omitting the independent variable of interest) for predictors of forewing length (mm). While forewing length tends increase with increasing developmental temperatures (non-significant correlation) as well as collection date, developmental temperatures and collection date interact. Forewing length increases with year One striking result from our studies is that mean forewing length, wing melanism, and setal length in C. meadii have increased, not decreased, during the past 60 years in this region. These phenotypic changes have occurred despite overall climate warming over this time period: mean July maximum temperature increased significantly from 1953 to 2013, (p < 0.001; slope estimate (se) = 0.03 (0.008)°C), resulting in a total increase of~1.8°C over 60 years (Additional file 1: Figure S1A). These historical shifts are inconsistent with simple expectations for evolutionary responses to climate warming. However, previous analyses suggested complex patterns of selection: increases in melanism are consistent with our expectations for high elevation populations that are flight limited [26]. Previous modeling of C. meadii also highlighted the importance of inter-annual variability [26]. Such environmental variation generates fluctuations in the magnitude and direction of selection that strongly reduce the rate of adaptive evolution to directional environmental change [26,46].
Several lines of past and current evidence support the contribution of phenotypic plasticity to variation in wing melanism. Lab experiments show that lower temperatures experienced during pupal development increase adult wing melanism in montane Colias [29,36,47]. Our historical analyses indicate that greater wing melanism is associated with lower pupal temperatures, and with earlier collection dates during the flight season (when temperatures are on average lower). These patterns are consistent with the idea that phenotypic plasticity in wing melanism in response to variation in pupal temperature is contributing to the historical pattern in this trait.
Whether such plasticity is adaptive depends on whether environmental cues during pupal development are good predictors of environmental conditions experienced by adults [48][49][50]. June and July temperatures in this region were only weakly correlated during the past 60 years (r = 0.19, p = 0.18), suggesting that pupal temperature in June is often a poor predictor of temperatures during the flight season during July (Additional file 1: Figure S2). When we consider the estimated pupal temperature for the specimens included in the analysis, there is a slight but significant increase in mean pupal temperature over the whole 60 year time period, but the regression explains only 3.2% of the variation in pupal temperature (slope = 0.015, R 2 = 0.032), with tremendous variability among years (residual variance = 5.50°C 2 ). This interannual variability may be contributing to the We present coefficients for predictor variables with corresponding standard errors, Z-scores, and the probability that the Z-score deviates significantly from 0. Note developmental temperature is a 40 day window encompassing larval and pupal development observed increase in melanism and setae length by causing variation in the direction of selection and mismatching developmental and flight temperatures. Thermal extremes that cause overheating are relatively rare in our montane sites [16], their incidence may be mostly independent of mean active season temperatures, and variation in thermoregulatory traits may be relatively ineffective in avoiding overheating. Thus, selection to increase flight time in cold years may be stronger and more consistent than selection to avoid overheating in warm years.
Change in the distribution of collection dates of historic specimens can indicate shifts in phenology [51]. We detect a phenological shift toward later collection in later years that interacts with elevation: individuals were collected earlier in the season at lower elevations and later in the season at high elevation in later years (Fig. 4). Climate warming may have extended the active season of high elevation individuals. The dependence on elevation may also reflect greater opportunity for phenological advancement at lower elevations where snow melts earlier. Exposure to cool, early or late season temperatures during development could select for wing darkening and setae lengthening and extend the developmental duration, increasing body size. Phenological shifts could result in the collection of individuals with darker wings and longer setae due to plastic responses to experiencing colder conditions earlier in the season at low elevation or cool late season conditions at higher elevations. The influences of phenological shifts at low and high elevation are difficult to distinguish.
The phenological shift toward later collection in later years at high elevation may indicate individuals from high elevations increasing fitness relative to lower elevation individuals as climate warms. These shifts could result in a higher proportion of individuals with high elevation phenotypes being sampled over time. Specimens in collections are concentrated over a  We present coefficients for predictor variables with corresponding standard errors, Z-scores, and the probability that the Z-score deviates significantly from 0 sufficiently narrow elevation range that we do not detect the elevational gradient in melanism that is well established [20,52]. The only significance effect of elevation is the interaction with year as a predictor of collection date. Fitness increases at higher elevation could have resulted in dispersal of individuals with high elevation phenotypes to lower elevations. Our results suggest that any evolutionary responses in Colias butterflies have been more complex than a steady shift toward lighter wings, shorter setae, and smaller body size in response to warming. Factors such as mate choice [53] or interspecific competition may also contribute to selection on wing coloration [54], further complicating temporal trends. Collection biases may have also contributed. Collected individuals are biased towards those in flight, but the importance of flight to fitness suggests that collections are consistently biased toward individuals with trait values associated with higher fitness.
Related studies suggest that our findings may be broadly applicable to understanding how thermoregulatory traits mediate interactions with the environment. An analysis of images in field guides for 473 European butterfly and dragonfly species indicated that cooler climates contain darker species and that distribution shifts between 1988 and 2006 resulted in assemblages becoming lighter in regions that warmed [11]. However, whether wing coloration has any thermoregulatory significance was unclear for most species in the analysis [55,56]. Although increasing melanism in response to decreasing developmental temperatures has been widely reported in insects, this includes many cases in which melanism plays no role in thermoregulation [57][58][59][60]. Our findings correspond to those of J Stamberger [34] who documented decreases in wing melanism and setal length in C. meadii with warmer temperatures, but no trend across years, for a more restricted time period and set of sites in Colorado.

Conclusions
Our finding that thermoregulatory traits and body size have increased over time highlights the complexity of biological responses to climate change. Environmental variability within and across seasons can limit or alter evolutionary responses. Variation in selection across spatial climatic gradients can also result in trait shifts, particularly if warming differentially impacts fitness across climatic gradients. Phenological shifts can also alter selection. Our results also highlight the importance of considering how plasticity can modulate species' responses to climate change. Temperature variability both within and across seasons and the potential for phenological, fitness, and distributional shifts may lessen the potential for responding to climate change via adaptive evolution.

Additional file
Additional file 1: Figure S1 Figure S2. The relationship between June mean temperatures and July maximum temperatures for each year from 1953 to 2013 at the weather station in Climax, CO (39.37 N, 106.18 W, [61]). Figure S3. The distribution of collection dates for the measured specimens. Table S1. Results Year: Elevation Partial residual of collection date (J) Fig. 4 Partial residual plot for normalized predictors of collection date (J). Collection date advances with year and elevation. However, the predictors interact such that collection date has shifted earlier in the season at lower elevation and later in the season at higher elevations in later years of averaging the top spatially autoregressive maximum likelihood models for wing length, wing melanism, and setae length. The top (by AICc) 20 subsets of the full model are averaged across the 50 iterations of the analysis We present coefficients for predictor variables with corresponding standard errors, Z-scores, and the probability that the Z-score deviates significantly from 0. Importance indicates the percent of top models that include each predictor variable. (DOCX 435 kb) Abbreviations T: Temperature in degrees Celcius; SAR: Spatially autoregressive; AIC: Akaike Information Criterion