Altitudinal trends in the phenology of butterflies in a mountainous area in central Spain

One of the best documented effects of climate change on biodiversity are shifts in phenology. However, long-term data quantifying and projecting the expected changes in phenology associated with climate warming are limited to a few well-recorded areas in the world. In the absence of temporal recording, an alternative approach is to determine the phenological response of species along marked gradients in climate or along latitudinal or altitudinal transects (space-for-time substitution). We studied the phenology (timing and duration of the flight period) of butterflies in 2006 along an altitudinal gradient (900–1680 m; estimated temperature lapse rate = –6.6°C/km) in the Serranía de Cuenca (central Spain) at the assemblage and individual species levels. Timing of the flight period was later for assemblages at high than at low altitudes. A similar trend of an increasing delay in the flight period with altitude was recorded for some individual species. However, there were also some exceptions to this pattern regardless of the number of sites and the altitudinal ranges of the species, suggesting possible local adaptation to regional climate. The duration of the flight period was shorter at high altitudes for assemblages, but this trend was not mirrored in the response of individual species. The results partly support substituting space-for-time when assessing the potential effect of climate change on phenophases such as the timing of the flight period, but we recommend extreme caution in extrapolating the results in the absence of information on how the responses of populations differ.


INTRODUCTION
Phenology has attracted increasing attention in recent years because it is particularly sensitive to climate change (Parmesan, 2006).Most phenological studies report marked temporal trends in numerous aspects of plant and animal life cycles associated with warmer temperatures (e.g., Walther et al., 2002).However, those surveys usually need to be carried out over long periods of time (several decades) if trends in phenology are to be detected (Parmesan, 2006).Consequently, there is a dearth of phenological information for many areas other than those intensively studied in Northern and Central Europe, and North America, particularly for a diverse and complex groups of organisms like insects (e.g., Sparks & Yates, 1997;Roy & Sparks, 2000;Forister & Shapiro, 2003;Hassall et al., 2007;van Strien et al., 2008;Westwood & Blair, 2010; but see Stefanescu et al., 2003;Gordo & Sanz, 2005, 2006 for southern European examples; and Doi & Takahashi, 2008 for an Asian example).Nevertheless, given the rapid rate in the increase of mean global temperature that has occurred during the last century and the projected rise for the coming decades (IPCC, 2007), there is an urgent need for research in other geographical areas, particularly those that are potentially more vulnerable to climate warming, like mountainous areas (Nogués-Bravo et al., 2007; see also Wilson et al., 2005).
An alternative approach to studying temporal variability in phenology is to study the climate determined changes in phenology that occur along altitudinal and lati-tudinal gradients, i.e., by using a space-for-time substitution (Hodkinson, 2005).Altitudinal gradients are particularly useful for this because they are usually characterised by marked changes in temperature (on average a 6°C decrease per km increase in altitude) over short geographical distances (a few tens of kms) (Körner, 2007).In addition, the effects of latitude-dependent parameters, like photoperiod, which also could potentially affect species phenology, are minimized (Fielding et al., 1999;Hodkinson, 2005).It is suggested that the responses of species to existing gradients in climate could be used to predict the likely effects of future climate change (Whittaker & Tribe, 1996;Fielding et al., 1999; but see Phillimore et al., 2010;Hodgson et al., 2011).
The few studies on the spatial (latitudinal and altitudinal) trends in insect phenology has revealed that poikilothermic organisms become active later in the colder areas at the highest latitudes and altitudes (e.g., Gutiérrez & Menéndez, 1998;Fielding et al., 1999;Roy & Asher, 2003;Merrill et al., 2008;Ashton et al., 2009;Gordo et al., 2010).This is to be expected because developmental rates are temperature-dependent (Bale et al., 2002).However, there are some significant exceptions, which show no clear trends in the timing of appearance, even for large data sets (e.g., Brakefield, 1987;Pollard, 1991;Roy & Asher, 2003; see also the early work by Verity, 1920).In addition, other important phenological attributes, such as the duration of the activity period (Memmott et al., 2007), are rarely assessed in spatial studies (but see Brakefield, 1987;Pollard, 1991;Gutiérrez & Menéndez, 1998;Kovanci & Kovanci, 2006).Therefore, further evaluations of the spatial trends in the phenology of insects (including both timing and duration) are needed for a better understanding of their potential susceptibility to climate change.In this sense, butterflies are a particularly appropriate group for studying phenological changes linked to environmental gradients, because, as poikilothermic animals, their life cycle is markedly affected by climatic conditions (Dennis, 1993).
In this paper, we document the phenology of butterflies at different altitudes in the Serranía de Cuenca (Central Spain).This is a species-rich mountain region in Spain, which is a prime butterfly area in Europe (Munguira et al., 2003), containing a relatively large number of Iberian endemics and endangered species (de Arce Crespo et al., 2004Crespo et al., , 2006)).We firstly studied the altitudinal patterns in timing and duration of the flight period of whole assemblages of butterflies and then that of the individual species.We expected a gradual delay in the timing and a progressive shortening of the duration of the flight period with increase in altitude.

System studied
Field work was carried out in the Serranía de Cuenca, Sistema Ibérico (province of Cuenca, central-eastern Spain).The Serranía de Cuenca is a mountain range of 7275 km 2 , running from 40°5´N 2°8´W (UTM grid reference 30TWK73) in the south-west to 40°20´N 1°43´W (UTM grid reference 30TXK06) in the north-east.This protected natural area ranges in altitude from c. 900 to 1866 m (Alonso Otero, 1998).The climate is typically continental Mediterranean, with wide daily and seasonal variation in temperature, and hot summers and cold winters.
Survey sites were open areas occurring in natural or seminatural habitats (usually woodland clearings, shrub or pasture) selected on the basis of accessibility and to provide a representative sample of the range in altitude in the region (Fig. 1; see Appendix for details).On the first visit to each site, a 500 m transect passing through the typical habitat of the location was established, and Universal Transverse Mercator (UTM) coordinates were recorded to the nearest metre at least every 100 m using a handheld Cobra GPS unit.Each transect was plotted in ArcView 3.2 (ESRI, 2001).
Butterflies were recorded at 10 sites in 2006 (Fig. 1).The mean distance between the sites that were closest to one another was 5.2 ± 0.7 (SE) km so for the species analysed transects represent independent populations in which phenological measures depend much more on local emergence than on immigration from neighbouring sites.Standardized 500 m long × 5 m wide transects were walked at each site every two weeks from May to September when conditions were suitable for butterfly activity (Pollard & Yates, 1993), On some occasions, stops were made to resolve problems over identification.A few butterflies were caught with a net and identified in the laboratory using field guides (Fernández-Rubio, 1991;Tolman & Lewington, 2002).Male genitalia were removed from those species that are more difficult to identify (e.g.Pyrgus spp.) and used to identify the species using Higgins (1975) and Fernández-Rubio (1976, 1977, 1981, 1982).Nomenclature follows García-Barros et al. (2004).

Butterfly phenology
The two most often used descriptors of activity patterns are timing and duration, or synchronisation (Wolda, 1988).Because temperate species show an approximately normal frequency distribution in the numerical occurrence of individuals during the course of a season, we described activity patterns using the mean flight date (hereafter mean date) and standard deviation about mean date, respectively (hereafter SD about the mean date) (Brakefield, 1987).Mean date is one of the most accurate estimates of shifts in phenology of organisms with unimodal temporal distributions (Moussus et al., 2010).Date was measured in terms of the number of days that elapsed since 31 March (1 = 1 April).Analyses were performed at two levels, the whole assemblage and individual species.At the level of an assemblage, we considered all the species and individuals recorded at each site.At the level of a species, because mean date and SD about the mean date are only appropriate for patterns with a clear activity peak, we analysed data only for those species with one annual generation (univoltine).In order to ensure that analyses included species whose phenology was representatively sampled, we considered only those species for which two or more individuals were recorded at least three transect locations, and excluded those sites where only one individual was recorded.We also excluded: (a) spring species for which the complete phenology was not recorded (e.g., Anthocharis cardamines), (b) species that overwintered as adults (e.g., Gonepteryx rhamni), and (c) species which presumably aestivate in summer in central Spain (e.g., Hipparchia semele) (García- Barros, 1987Barros, , 1988)).Information on voltinism was based on published records (Fernández-Rubio, 1991;Tolman & Lewington, 2002) and more specific studies in the study area (de Arce Crespo et al., 2004, 2006, 2009).
Changes in timing and duration of the flight period at the assemblage and species levels were tested by regressing mean date and SD about the mean date, respectively, against the altitude of each transect site.To determine whether the slopes of regressions tended to show the same sign for a significant percentage of species, we used two-tailed sign tests.
Because of unsuitable (cloudy) weather conditions, some transects could not be walked during some of the two-week  1. Results of the linear regressions of mean date (days elapsed since 31 March, 1 = 1 April) against altitude (km) for 33 univoltine species.The species are arranged in alphabetical order and those with significant regressions are in bold.The number of sites occupied, the minimum and maximum altitudes (m), and the total number of individuals recorded (n) for each species are also shown.
periods: localities 1, 3, 5, 7 and 9 on 13 September, and locality number 1 on 27 September (Appendix).In those cases, we estimated species abundances for the missing periods by calculating the arithmetic mean between the counts in adjacent periods following Pollard & Yates (1993).To control for the potential biases that resulted from missing sampling periods, all analyses were performed on the raw data (including transect gaps) and then by using the data that included the estimated abundances of the species.

Lapse rate
We determined the change in temperature along the altitudinal gradient in the study area (lapse rate).We selected all those meteorological stations in the area with complete records of temperature, including records for the year studied: 6 stations in Cuenca, Guadalajara and Teruel provinces ranging between 956 and 1601 m in altitude in 2001 and 2004 to 2006.Lapse rate was calculated by regressing mean annual temperature on altitude (km).
All the statistical analyses were performed using SPSS software (SPSS Inc., 2001).

RESULTS
A total of 111 species and 7555 individuals was recorded at the 10 localities visited.Thirty-three species satisfied the criteria for analysis (Table 1).

Phenology of assemblages
There was a 19 day-difference between the mean dates for the butterfly assemblages at lowest and the highest sites.This delay was reflected in a positive and significant relationship between the mean date for assemblages and altitude [R 2 = 0.60, P = 0.009, n = 10; mean date (days since 31 March) = 80.64 + 24.48*elevation (km)], which indicates a delay of c. 24 days in the phenology of butterflies for every km increase in altitude (Fig. 2).Using data that included estimates of the abundances for missing sampling periods led to similar results, with a slightly longer delay of c. 28 days per km [R 2 = 0.634, P = 0.006, n = 10; mean date = 77.21+ 27.61*altitude (km)].
The relationship between the SD of the mean date and altitude was negative and significant for assemblages [R 2 = 0.433, P = 0.039, n = 10; SD of the mean date (days) = 40.10-8.80*altitude (km)].This indicates a progressive shortening of the butterfly flight season with increase in altitude of c. 9 days per km (Fig. 2).Using data with esti-mated abundances increased the scatter of the relationship (but the magnitude of change was similar), but it was not significant at P < 0.05 [R 2 = 0.321, P = 0.088, n = 10; SD about the mean date = 42.84-9.62*altitude (km)].

Phenology of species
For the 33 univoltine species the slopes of the regression between mean date and altitude for 23 were positive and for 10 they were negative (mean ± SE, 18.50 ± 7.84; sign test: P = 0.037, n = 33) (Table 1).Six regressions were significantly positive (A.crataegi, C. arcania, A. adippe, M. lachesis, P. tithonus and T. sylvestris), indicating a delay in the flight period with increase in altitude and none of the negative regressions were significant.The magnitude of the delay ranged from 16 days per km for M. lachesis to 54 days per km for T. sylvestris (mean ± SE, 39.13 ± 5.91).Using data that included estimated abundances led to similar results, with 24 positive slopes and 9 negative slopes (sign test: P = 0.015, n = 33) (only the non-significant regression slope for H. statilinus changed from negative to positive; results not shown).Apart from the same 6 species cited above, Pyrgus cirsii also had a significantly positive regression between mean date and altitude.
Table 2 shows the results for the regressions between the SD about the mean date and altitude.Fifteen species had positive slopes and 17 negative slopes (no regression parameters were estimated for C. glycerion due to zero variance in the data) (mean ± SE, -1.58 ± 3.85; sign test: P = 0.860, n = 32).One regression was significantly positive (E.triaria) and one significantly negative (A. arethusa) (mean ± SE,29.54 ± 55.22).Using data that included estimated abundances produced 16 positive and 16 negative slopes (excluding also C. glycerion) (sign test: P = 1, n = 32).The regression for E. triaria remained significant but not that for A. arethusa.In contrast, the regression for C. briseis became significant when estimated data were included in the analysis (results not shown).

Timing of appearance
The temperature lapse rate recorded for the altitudinal gradient sampled (-6.6°C/km) was similar to that recorded in other mountainous areas in central Spain (e.g., Wilson et al., 2005) and elsewhere (Hodkinson, 2005).At the assemblage level, there was a c. 24 day delay in the mean date per km increase in altitude (Fig. 2).Thus, butterfly assemblages emerged sooner at low altitude and warm sites than at cold high altitude sites, as recorded using a similar approach in other areas (e.g., Gutiérrez & Menéndez, 1998).For individual, univoltine species, there was also a significant trend for positive slopes for the altitudinal delay in mean date, and all significant slopes were positive, suggesting a similar pattern at the species level.Nevertheless, except for M. lachesis (c.17 days/km), all significant slopes for individual species (c.31-54 days/km) were steeper than those for the whole assemblages.The altitudinal delay for A. crataegi in the Serranía de Cuenca (33.08 days/km) was extremely close to that recorded in the Sierra de Guadarrama (also in central Spain) in the same year using an identical approach (33.1 days/km, Merrill et al., 2008), which supports the robustness of our results.Nevertheless, it is very likely that altitudinal delays in mean date change from year to year depending on conditions (e.g., Merrill et al., 2008;Gutiérrez Illán et al., unpubl. results) and these rates should not be taken as fixed attributes for these species.
The altitudinal delay at the assemblage level is very likely the result of the sum of the individual altitudinal delays of a significant proportion of the species and also of the different contributions of species that emerge late in a season along the altidinal gradient.Two of the most abundant species in the area, E. zapateri and P. coridon (1376 individuals in total, Table 1), emerged in late summer and occurred mainly at high altitudes, which could have resulted in the delay in mean flight dates of whole assemblages recorded at high altitudes.
Delays in the time of appearance or other phenophases (such as flowering, leaf unfolding or autumn colouring) associated with altitude, although expected, are rarely reported (e.g., Gutiérrez & Menéndez, 1998;Fielding et al., 1999;Blionis et al., 2001;Schleip et al., 2009;Gordo et al., 2010).In the case of insects, delays in their appearance at high altitudes may be due to the direct effects of temperature on the time taken for eggs to hatch (e.g., Whittaker & Tribe, 1996;Fielding et al., 1999) and on the developmental rate of immature stages (e.g., Bale et al., 2002), or to indirect effects of resource quality on growth (e.g., Erelli et al., 1998;Suzuki, 1998).Another possibility is that timing of appearance is correlated with the termination of winter diapause, which can be controlled by factors correlated with altitude, such as temperature and moisture (Leather et al., 1993).
However, no evidence at all was found for a delay in timing of appearance associated with altitude for at least 10 (30%) species (those with negative slopes -the most conservative estimate, Table 1), which contrasts with pre- vious studies on altitudinal gradients (Gutiérrez & Menéndez, 1998 report only 1 out of 41 regressions have a negative slope).Interspecific variability in the magnitude of the delays did not appear to be associated with a large sample or a wide range in altitudes, because regressions between altitudinal delay and number of sites and altitudinal delay and altitudinal range were nonsignificant (P = 0.916 and P = 0.675, respectively).It is suggested that spring events are more affected by climate than those occurring later in a season (Menzel, 2000;Forister & Shapiro, 2003).This point has been poorly studied along spatial gradients (but was mentioned by Verity, 1920), but we have found no evidence for this in our data based on regressing the average mean date on altitude (average for all sites where one particular species was present) (P = 0.256).Nevertheless, our study does not include the species that emerge very early in a season because the study started mid-May, which could have had a strong effect on the relationship between average mean date of emergence and altitude.Synchrony in mean flight date across temperature gradients tends to indicate adaptations of butterfly populations to local regional climates (Roy & Asher, 2003).Such adaptations could include behavioural and physiological processes, like selection of warm microhabitats for egg-laying at high altitudes, or life history adaptations in terms of larval growth rate or adult body size (Hodkinson, 2005).

Duration of the period of activity
The results show a decline in the duration of the flight period of butterfly assemblages with increase in altitude (Fig. 2), which accords with previous results reported for northern Spain (Gutiérrez & Menéndez, 1998).However, we found no evidence for a similar pattern in the flight periods of individual species (Table 2): linear regressions of the SD about the mean date against altitude resulted in similar proportions of positive and negative slopes, considering only the significant relationships.Thus, the pattern at the assemblage level was not the result of the simultaneous, progressive shortening of the flight period of individual univoltine species.Shortening of the flight period of the whole assemblage with increase in altitude (with no apparent shortening of the flight periods of individual species) could arise if species flying later in the season showed little response to increase in altitude, but, as stated above, we found no evidence for this.An alternative explanation is that the multivoltine species complete fewer generations per year at the higher altitudes, but there was no evidence for this either because the progressive shortening with increase in altitude was also detected at the assemblage level when only the 33 univoltine species were included in the analysis (P = 0.032).

Relationship of phenological changes to climate change
Phenological responses of species along an altitudinal gradient can be used to predict the likely effects of possible future changes in climate (Whittaker & Tribe, 1996;Fielding et al., 1999).In 2006, the mean date was delayed c. 24 days per km for butterfly assemblages and c. 17-54 days per km for univoltine species for which the relationships were significant (Table 1).Over the period 2001-2006, the annual mean temperature decreased by c. 6.6°C per km increase in altitude (Fig. 3).This suggests that a 1°C rise in mean temperature could advance butterfly phenology by 3.7 days for entire assemblages and 2.6-8.2days for individual univoltine species in this region.These trends are broadly in line with those found by Roy & Sparks (2000), who report an advance of 2-10 days/°C for the mean date of appearance of British butterflies.
Based on different scenarios proposed for the late twenty-first century, climate models project global temperature increases ranging 1.8-4.0°C(IPCC, 2007).This would advance butterfly phenology by 6.7-14.8days for whole assemblages and by 4.6-32.8days for individual univoltine species.If the appearance of interacting species (host plants and natural enemies) shifts at different rates, this would lead to a mismatch in phenology that may have severe consequences for butterfly population dynamics (Visser & Both, 2005).
Nevertheless, there was a large proportion of species for which the mean flight date is synchronized along the altitudinal gradient.This raises doubts about the validity of using the space-for-time substitution approach in studies in which spatial replication is possible but temporal replication is much more difficult (Phillimore et al., 2010).In the absence of knowledge on the temporal relationship between phenology and temperature, estimates for projected shifts in appearance based on spatial data should be made with extreme caution.Otherwise, the magnitude and even the sign of rates of phenological change with climate could be equivocally estimated (Phillimore et al., 2010).In this sense, further research testing the consistency between rates of phenological shift based on temporal data vs. those based on spatial data will be particularly useful.
Our work suggests that, for whole butterfly assemblages and some individual species, there is a delay in phenology with increase in altitude, which is supportive of using climatic gradients to evaluate the likely effect of climate change on phenological shifts.However, because a substantial proportion of species appeared to emerge synchronously along the altitudinal gradient, the spacefor-time substitution approach in phenological studies should be used with extreme caution.In cases like this, some evidence that there is a consistent relationship between spatially-and temporally-based rates of phenological change is desirable.

Fig. 1 .
Fig. 1.Map of the study area showing the altitude and locations sampled (white triangles).200 m altitude bands are shown with increasing shading from < 900 m to > 1700 m a.s.l.The border of the protected environmental area in the Serranía de Cuenca is indicated by a white line.The numbers next to the triangles correspond to those in Appendix.

Fig. 3 .
Fig. 3. Relationship of mean annual temperature (°C) with altitude (km) for the period 2001, 2003-2006 based on temperature records of meteorological stations located in the Sistema Ibérico.

TABLE 2 .
). Results of the linear regressions of SD about the mean date against altitude (km) for 33 univoltine species.The species are arranged in alphabetical order and those with significant regressions are in bold.No regression parameters were estimated for C. glycerion due to zero variance in the data.