Because they are small yet easily observed, butterflies have been the subject of numerous studies of the effects of climate change on living organisms (Dennis 1993; Parmesan 2006). The responses have included earlier emergence (Roy & Sparks 2000; Forister & Shapiro 2003; Stefanescu et al. 2003; Kearney et al. 2010; Diamond et al. 2011), shifting abundance and range (Parmesan et al. 1999; Hellman et al. 2008; Breed et al. 2013), expansion of larval diet (Pateman et al. 2012), and extirpation of populations (McLaughlin et al. 2002). The many influences of climate change on butterflies, as well as on other organisms, are still being discovered.

A study by Polgar et al. (2013) assessed the influence of temperature on the appearance of 10 univoltine lycaenid butterfly species, including five spring emergent species (Callophrys spp., elfins) and five summer emergent species (Satyrium spp., hairstreaks). Data for their analyses came from the records of the Massachusetts Butterfly Atlas and the Massachusetts Butterfly Club, an active citizen science organization that has compiled extensive observation records from the 1990s on. Citizen science data like these are becoming more important (Dickinson et al. 2010) because they provide abundant data over broad areas and, as a result, are being used increasingly in scientific studies. Because of the extensive records of butterflies accumulated locally by citizen scientists, Massachusetts is the U.S. state with the most records available for studies like those of Breed et al. (2013), Polgar et al. (2013), and the one presented here.

This paper extends the analysis of Polgar et al. (2013) of 10 species of lycaenids (listed in Table 1) in four ways: we have (1) expanded the data set by adding observations from recent years; (2) applied a different analytical method (Hassall & Thompson 2010) that minimizes possible bias from unequal sampling over time; (3) compared the analysis of observational data with and without equal sampling; and (4) calculated more statistically significant estimates of the rates of phenological advance for these species. We expected that, in response to climate warming, all 10 lycaenid species would appear significantly earlier in recent years than they did in the 1980s; that the phenological advancement of these species would be greater than that estimated by Polgar et al. (2013) because of the warmth of recent years; that the use of equal sampling for analysis would provide more reliable estimates of phenological change; and that spring emergent species would show greater phenological advancement than summer emergent species.

Methods

We began with the set of citizen science records from 1986 to 2009 analyzed by Polgar et al. (2013), to which we added data for years 2010 through 2012 obtained from the Massachusetts Butterfly Club ( http://www.massbutterflies.org), as reviewed by one of us (SBS). After processing the data in different ways (described below), we analyzed phenological shifts by first regressing observation dates against year for each species separately. Then for the groups of elfins (early) and hairstreaks (late), we regressed observation date against year, against average spring (Mar-May) temperatures, against summer (Jun-Aug) temperatures, and simultaneously by multiple regression against both year and spring temperature. To calculate the days of phenological advancement for each species or group, we multiplied the slopes of the regressions (days per year) by 26 because 2012 followed 1986 by 26 years.

From the full set of observation records over the 27 year period, we chose two sets of data for the regression analyses. The first included all observations, while the second followed Polgar et al. (2013) in using only the first 20% (FTP) of observations within each year (with a five observation minimum). The FTP procedure captures the clearest responses to climatic conditions by focusing on the timing of the initial emergence of each species without the extended tail of later emergence; it also eliminates any bias that may develop from observers paying less attention to a species later in its flight period.

In the analysis of these two data sets (All and FTP), we followed the recommendation of Hassall & Thompson (2010) to eliminate possible bias from uneven recorder effort over time by analyzing the same number of records in different time periods. Results may be skewed by unequal sampling effort, and it is often the case that collection records are more numerous in recent years. Therefore, we divided 1986–2012 into nine consecutive three-year periods, and for each species we found the smallest number of observations in any of the three-year periods. Then, using that smallest number to be the usable number of records available for each three-year period (with a three record minimum), we randomly resampled data from the more frequently recorded time periods so that equal numbers of records were analyzed for each threeyear period. Random sampling was accomplished by assigning spreadsheet-generated random numbers to each observation record, sorting those records numerically by random number within each three-year period, and choosing the first usable number of records for each period. The data sets with equal sampling are labeled All-equal when all observations were used and FTP-equal when FTP data were used. Using both data sets, we ran three randomizations separately for each species as well as for elfins and hairstreaks as groups and then averaged the results.

Climate change in Massachusetts from 1986 to 2012 was assessed by examining the mean monthly temperatures for three different weather stations using data obtained from the National Climatic Data Center ( www.ncdc.noaa.gov). The sites chosen for analysis represent different regions of the state: Plymouth lies in the lowland of eastern Massachusetts; Worcester is on the Worcester plateau near the middle of the state; and Amherst is in the Connecticut River valley towards the west. Combined, these three sites give a measure of average temperature across the state and represent regions with many of the butterfly observations. We analyzed mean monthly temperatures for March, April, and May (spring months) and June, July, and August (summer months). To assess changes in climate, we regressed mean temperatures for each month and for all six months combined from these three sites against year. All statistical analyses were accomplished with SPSS Version 21 (IBM Corp. 2012), with Zar's (1996) method to compare slopes of different regressions.

Results

Earlier emergence was apparent and significant in all 10 lycaenid species we examined (Table 1, Fig. 1). Furthermore, advancement was greater for springemerging elfins than for summer-emerging hairstreaks (Figs. 1 & 2). From 1986 to 2012, based on FTP-equal data, advancement of the flight season for elfins ranged from 10.0 days in Henry's Elfin (C. henrici) to 18.2 days in Eastern Pine Elfin (C. niphon), while the advancement for hairstreaks ranged from 5.3 days in Banded Hairstreak (S. calanus) to 9.6 days in Edward's Hairstreak (S. edwardsii). With all species combined into the elfin and hairstreak groups, emergence was 14.2 days earlier for spring-emergent elfins and 7.9 days earlier for summer-emergent hairstreaks, and the difference in slopes was highly significant (Fig. 2; t₄₄₂>100, p < 0.001).

Table 1.

Regression statistics for analysis of phenological advancement of 10 lycaenid species. Drawing from both All and FTP (first 20%) data, equal sampling was used from nine three-year time periods from 1986–2012. The average of slope (beta), std. error, r², and sample size from three randomizations is presented here, along with the p values from all three randomizations.

Table 2.

Comparison of regressions for all elfins and all hairstreaks against year from different methods for analyzing the data. All includes the set of all observations from 1986 to 2012; All-equal is based on the full set but with observations selected randomly to give an equal number of observations from each three-year period; FTP includes the first 20% of all observations from each year; and FTP-equal was based on only the first 20% of observations from which an equal number of data points were chosen for each three-year time period. A single randomization is shown for the equal analyses. For each regression, the slope (beta), r², sample size, and significance of the regression are shown.

The different methods of analysis yielded similar results, with all regressions producing significantly negative slopes that indicate earlier emergence (Tables 1 & 2), even with p=0.01 as the threshold for significance because of the use of multiple tests. Despite having smaller sample sizes, the regressions for each species using FTP-equal data were almost always significant, failing only with small data sets: Hoary Elfin (C. polios) (insufficient data for FTP-equal analysis); Frosted Elfin (C. irus) with one of three randomized analyses (23 data points); Henry's Elfin (C. henrici) with two of three randomized analyses (24 data points); and Acadian Hairstreak (S. acadica) with two of three randomized analyses (27 data points; marginal significance). Analyzing the larger samples from the Allequal data sets, similar results were obtained, again lacking significance with only some randomizations for Hoary Elfin (C. polios), Frosted Elfin (C. irus), and Henry's Elfin (C. henrici) (Table 1). Marginal significance also occurred with one randomization for Banded Hairstreak (S. calanus).

Using data combined for elfins and hairstreaks, we compared the two different methods of analysis: All vs. FTP, and equal vs. non-equal sample sizes (Table 2). It is apparent that the slopes of the regressions were very similar whether the analysis was of All data or FTP data and whether we chose equal sampling across three-year time periods for analysis or not. The sample sizes were smaller when using FTP data, but the regressions were stronger (larger r² values, with 2.4 to 4 times more variation accounted for by the regression). Equal sampling in our data set gave a steeper slope of the regressions for elfins but not for hairstreaks.

When the day of observation was regressed against average spring temperature, an even stronger result was obtained than from the regression against year (Fig. 3). Using FTP-equal data, phenological advancement was 4.80 days/°C for elfins and 3.14 days/°C for hairstreaks, with r² values of 0.305 for elfins and 0.290 for hairstreaks. In contrast, when regressed against year, the same two groups showed smaller r² values of 0.194 and 0.147, respectively. As expected, elfins responded more strongly to warmer conditions than did hairstreaks, and the difference in slope was significant (t₄₄₂ = 27.65, p < 0.001). Multiple regression provided a simultaneous evaluation of the relative strength of the effects of temperature and year. Both predictors were significant at p<0.001, with temperature having a greater effect than year. With FTP-equal data, the standardized coefficients for spring temperature were -0.451 and - 0.459 for elfins and hairstreaks respectively, while for year they were less at -0.275 and -0.213. Using multiple regression to compare the influence of spring temperatures with those of summer, we found that average summer temperature (from the months of June, July, and August) did not affect spring-appearing elfins (t₁₇₀=1.264, p=0.206), and though the effect of summer temperature was significant on hairstreak appearance, it was less than the effect of spring temperature (standardized coefficients, spring = -2.910, summer = - 1.462).

Fig. 1.

Phenological advancement of ten species of univoltine lycaenids from 1986 to 2012. The number of days of earlier emergence, shown on the y-axis, was determined using two different analyses of the observational data. All-equal was based on all observations; FTP-equal was based on the first 20% of observations within a year; and from these data sets, equal numbers of observations were randomly selected from each three-year period. The results of three such randomizations were averaged to produce the results given here (mean +/- SE). See text for further explanation.

Fig. 2.

Phenological advancement across years of 4 spring-emergent Callophrys spp. (elfins, shown in open circles; data were insufficient to include Hoary Elfin, C. polios) compared to that of 5 summer-emergent Satyrium spp. (hairstreaks, shown in solid diamonds), using first 20% of observations for all species in each year, with equal sampling of the data (FTP-equal data). Average advancement from 1986 to 2012 was 14.2 days for the 5 elfins and 7.9 days for the 5 hairstreaks. Statistics are given in Table 2. The slopes are significantly different (t₄₄₂>100, p < 0.001).

Fig. 3.

Appearance date of the same 4 elfin species (open circles; data were insufficient to include Hoary Elfin, C. polios) and 5 hairstreaks (solid diamonds) regressed against each year's average spring temperature (°C, averaged from mean monthly temperatures for March, April, and May for Amherst, Plymouth, and Worcester, MA). Based on first 20% data with equal sampling (FTPequal data). For elfins, slope = -4.801, r² = 0.305; for hairstreaks, slope = -3.141, r² = 0.290. The slopes are significantly different (t₄₄₂=27.65, p < 0.001).

Our analysis of temperature records confirmed that the climate in Massachusetts has warmed significantly from 1986 to 2012. Over those 27 years, the average increase for all six warm-season months across the three stations we examined was +0.74°C, but the increase varied by month. The largest and most significant increase was +2.21°C for the month of April (F_1,77=13.968, p<0.001), but we also found significant increases of +1.13°C for July (F_1,77=5.122, p=0.026) and +1.47°C for August (F_1,77=11.098, p=0.001). March (+0.99°C), May (+0.40°C), and June (+0.46°C) showed non-significant increases. Temperature varied more during the spring months (March-May, 8.26° +/- 4.96°C, mean + std. dev.) than during the summer months (June-August, 20.69°C +/- 1.78°C). The last year of this study, 2012, was the warmest ever recorded for the U.S. as a whole (NCDC 2013), but based on the mean monthly average temperatures for our three representative weather stations across the months of March, April, and May, the warmest springs were in order 2010 (10.92°C), 2012 (10.40°C), 1991 (10.31°C), and 1998 (9.44°C). These years provide the data to the farthest right in Fig. 3.

Discussion

The results reveal the strong response of these 10 univoltine species to a warming climate; the emergence dates of all species have advanced significantly during the relatively short period of 26 years. No matter how one looks at the data, analyzed with All or FTP data and with or without equal sampling, the signal of phenological change is apparent. Earlier emergence correlates directly with the documented increasing average temperatures in the state. Our analyses show that, as Miller-Rushing & Primack (2008) found for flowering plants, spring temperature is a stronger predictor of emergence than year because the effect of year is through increasing warmth over time. Overall, we found that the five elfins are emerging 4.8 days/°C earlier, whereas the five hairstreaks are emerging 3.1 days/°C earlier. These results are greater than the 2.4 days/°C found by Kharouba et al. (2013) but similar to the response for the same species of 5.5 days/°C for elfins (based on Mar-Apr temperatures) and 2.9 days/°C for hairstreaks (based on May-Jun temperatures) reported by Polgar et al. (2013).

While our results of elfins and hairstreaks were similar to those of Polgar et al. (2013) in the number of days advanced per °C, the warmth of recent years has greatly accelerated emergence dates in those years. The previous study found that elfins were emerging 7.6 days earlier over 24 years, while we found them emerging 14.2 days earlier over the 27 years leading to 2012. For hairstreaks, Polgar et al. (2013) calculated 3.2 days advance (over 24 years), and we found 7.9 years advance (over 27 years). The contrast shows a rapid response by these butterflies to warmer conditions, a strong response similar to that of flowering plants, which have also advanced greatly in phenology in recent warm years (Ellwood et al. 2013).

The different ways we used to analyze observational data gave similar but slightly different overall results. The use of the first 20% (FTP) of observations focuses on initial appearance while excluding the extended tail of emergence; the use of equal sampling avoids potential bias from data being heavily weighted to one time period. The cost of each of these methods is a reduction in sample size. When data are sufficient, however, use of both methods can be valuable, as was the case with our study of univoltine butterflies. Our sampling of equal numbers of records from nine three-year periods worked well and gave us sufficient data for detailed analysis. We found greater explanatory power (higher r² values) when using FTP data and, for elfins, a stronger regression (steeper slope) when using equal sampling. These results highlight differences that resulted from the choices of FTP data and equal sampling.

Size of the data set influences the results, and the phenological shift was statistically weakest with Hoary Elfin (C. polios), the species with the fewest records. The most variable results between methods occurred with Henry's Elfin (C. henrici) and Frosted Elfin (C. irus), species for which the number of observations were small, particularly for FTP-equal data (Table 1). Despite the variable results for these species, their phenology has advanced. The least phenological shift, on the other hand, occurred in Banded Hairstreak. Among the hairstreaks, this species uses walnuts (Juglans spp.) and hickories (Carya spp.) in addition to oaks (Quercus spp.) as host plants (Opler & Krizek 1984). The host plants of the other hairstreaks include willows (Salix spp.) and members of the Rosaceae (e.g., Prunus spp.) and Ericaceae (e.g., Vaccinium spp.), in addition to Quercus spp. (for Edward's Hairstreak, S. edwardsii) (Opler & Krizek 1984). Walnuts and hickories are among the last trees to leaf out each year (Lechowicz 1984), with walnuts showing sensitivity to winter chilling as well as spring warming (Pope et al. 2013). Altermatt (2010) has described the influence of larval diet on phenological sensitivity, and if Banded Hairstreak (S. calanus) is tracking host plant phenology, the difference in host plants might be a factor in its reduced sensitivity to rising temperatures.

Spring species are more sensitive to temperature changes because of the faster rate of warming in spring than in summer and the greater variance in temperature during spring months. In fact, the greatest increase in temperature we detected in our Massachusetts data was for the month of April. Greater sensitivity of spring emergent butterflies has also been documented by Forister & Shapiro (2003) in California, Diamond et al. (2011) in the U.K., and Kharouba et al. (2013) in Canada. The median appearance date of the elfin species, based on all sightings across the full 1986–2012 period, ranged from day-of-year 124 for Brown Elfin (C. augustinus) through 138 for Frosted Elfin (C. irus), whereas the median date of appearance for the summer-emergent hairstreaks ranged from 191 for Banded Hairstreak (S. calanus) to 195 for Coral Hairstreak (S. titus). The five elfin species studied here overwinter as pupae, ready for eclosion, whereas the five hairstreaks overwinter as eggs (Opler & Krizek 1984). Diamond et al. (2011) emphasized that species overwintering in more advanced stages are those that emerge earlier in the year. Species that overwinter as pupae must find plants for oviposition during their spring flight period, while those that overwinter as eggs must find plants for larval feeding. Thus, the sensitivity of spring and summer emergent species to warming may depend in part on responses of their host plants to changing climatic conditions.

The responses of butterflies to climate change can indicate how other groups of organisms may respond because, in many regions of the globe, butterflies are sensitive indicators of habitat change (Merckx et al. 2013). Hodgson et al. (2011) have shown that analyses like the one done here can be used to predict effects of warming for several years in the future. In addition to phenology, range is sensitive to warming; using data from the same geographic state, Massachusetts, Breed et al. (2013) documented shifts in range of numerous butterfly species toward more northern areas. The results of our study add to understanding the effects of climate change on the seasonal appearance of univoltine butterflies, present a comparison of different ways to analyze observational data, and provide measurements of the greater sensitivity of spring emergents to warmer conditions.