Butterfly phenology in Mediterranean mountains using space‐for‐time substitution

Abstract Inferring species' responses to climate change in the absence of long‐term time series data is a challenge, but can be achieved by substituting space for time. For example, thermal elevational gradients represent suitable proxies to study phenological responses to warming. We used butterfly data from two Mediterranean mountain areas to test whether mean dates of appearance of communities and individual species show a delay with increasing altitude, and an accompanying shortening in the duration of flight periods. We found a 14‐day delay in the mean date of appearance per kilometer increase in altitude for butterfly communities overall, and an average 23‐day shift for 26 selected species, alongside average summer temperature lapse rates of 3°C per km. At higher elevations, there was a shortening of the flight period for the community of 3 days/km, with an 8.8‐day average decline per km for individual species. Rates of phenological delay differed significantly between the two mountain ranges, although this did not seem to result from the respective temperature lapse rates. These results suggest that climate warming could lead to advanced and lengthened flight periods for Mediterranean mountain butterfly communities. However, although multivoltine species showed the expected response of delayed and shortened flight periods at higher elevations, univoltine species showed more pronounced delays in terms of species appearance. Hence, while projections of overall community responses to climate change may benefit from space‐for‐time substitutions, understanding species‐specific responses to local features of habitat and climate may be needed to accurately predict the effects of climate change on phenology.

The most robust forecasts of phenological responses to climate change combine high-quality monitoring data collected over an extended time period (de Arce Crespo & Gutiérrez, 2011;Banet & Trexler, 2013) with an identified dependence of a phenological trait on reliable environmental cues (Reed, Waples, Schindler, Hard, & Kinnison, 2010). Detailed phenological time series have been analyzed for butterflies in northern and central Europe (Altermatt, 2012;Roy & Sparks, 2000;Van Strien, Plantenga, Soldaat, Swaay, & WallisDeVries, 2008) and elsewhere (Brooks et al., 2017;Diamond et al., 2011), but many other regions lack long enough time series of When no long-term data are available, an alternative approach to studying phenology is to substitute space for time, by assuming that the spatial relationship between an environmental factor (e.g., elevation) and a phenological response (e.g., time of appearance) can be used as a proxy for the temporal relationship (Banet & Trexler, 2013).
In this way, studies investigating how phenotypic traits change along latitudinal or elevational gradients can contribute to the prediction of species responses to climate change (de Arce Crespo & Gutiérrez, 2011;Gutiérrez & Menéndez, 1998;Hodkinson, 2005;Leingärtner, Krauss, & Steffan-Dewenter, 2014;Merrill et al., 2008). However, space-for-time substitution can become less valid at certain spatiotemporal scales (Blois, Williams, Fitzpatrick, Jackson, & Ferrier, 2013) or lead to underestimations of changes in diversity (França et al., 2016) especially under the pressure of a changing environment (Damgaard, 2019). Therefore, the implicit use of space-for-time substitution should be treated with caution in modeling community responses to climate change.
Information on species' ecological and life-history traits has also been used to further understanding of interspecific variation in phenological responses to climate change (Diamond et al., 2011;Kharouba, Paquette, Kerr, & Vellend, 2014;Leingärtner et al., 2014). For example, species whose flight periods occur earlier in the year and less mobile species appear more sensitive to temperature variation than late flying or more mobile species (Kharouba et al., 2014), and multivoltine species may be more able to respond to warming by increasing the frequency of their annual generations (Altermatt, 2009(Altermatt, , 2010. Although earlier emergence dates and increased numbers of generations have been widely documented, it has also been shown that some insects could be negatively affected by warmer climates. If juvenile stages complete development in late summer instead of entering the overwintering stage, a lack of sufficient time and suitable conditions to breed could lead to population declines ("the lost generation hypothesis") (Glazaczow, Orwin, & Bogdziewicz, 2016;Van Dyck, Bonte, Puls, Gotthard, & Maes, 2015;van der Kolk, WallisDeVries, & Vliet, 2016). Also for butterflies performing a photoperiodically induced summer dormancy, like Mediterranean Maniola butterflies (Van Dyck et al., 2015), climate warming might have negative effects on populations. If the summer drought became extended and the butterfly deposited her eggs before the onset of vegetation regrowth, triggered by a shortened photoperiod at the beginning of autumn, the young larvae would have no suitable fresh grasses to feed on and starve to death. Other negative consequences of climate change could include phenological mismatches between trophically interacting species such as butterflies and their host plants (Bale et al., 2002;Parmesan & Yohe, 2003;Visser & Both, 2005). Overall, climate change is responsible for well-studied phenological shifts, but their magnitude and direction can largely vary even between species inhabiting the same latitude (Diez et al., 2012). For these reasons, it is important to assess the effect of climate change on phenology from a comprehensive range of environments and ecological communities.
Elevational gradients are potentially useful space-for-time proxies because they combine significant variation in temperature over short geographic distances (Körner, 2007) with minimal variability in photoperiod (Fielding, Whittaker, Butterfield, & Coulson, 1999;Hodkinson, 2005). In addition, microclimate and habitat conditions (including vegetation structure and canopy cover) vary over elevational gradients (Suggitt et al., 2011) and can buffer ecological communities against coarse-scale trends and patterns in climate change (Gillingham, Huntley, Kunin, & Thomas, 2012;Suggitt et al., 2012). Therefore, phenology can vary markedly over elevational gradients but also within an altitudinal belt depending on habitat type, and Altermatt (2012) showed that the seasonal appearance of butterflies is influenced by both of these variables.
Testing local effects of elevation and habitat on phenology using space-for-time assumption could be a valid approach to understanding and predicting ecological responses to climate change (Banet & Trexler, 2013;Hodgson et al., 2011;Leingärtner et al., 2014). Although this method has some caveats (e.g., it cannot track year-to-year changes in species phenology), it can, however, serve as a short-term "tracking device" that mimics the longer seasons and milder winters that are expected as the climate warms (EEA, 2017;van der Wiel, Kapnick, & Vecchi, 2017).
In general, increasing elevation is expected to influence species' phenology by shortening the annual activity window, forcing stages in the life cycle to appear later while maintaining synchrony with resources and suitable environmental conditions (Brown & Lomolino, 1998;Despland, Humire, & Martín, 2012;Hodkinson, 2005). There is much evidence that temperate butterflies become active later annually at higher elevations (de Arce Crespo & Gutiérrez, 2011;Illán, Gutiérrez, Díez, & Wilson, 2012;Merrill et al., 2008;Shapiro, 1975). In addition, there is evidence that climate warming has led to both earlier appearance and extension of the flight period at high elevations (Konvička, Beneš, Čížek, Kuras, & Klečková, 2016). In this paper, using space-for-time inferences we frame phenological responses to climate change in a biota which lacks long-term data. By examining the phenology of butterfly communities in two mountainous areas in Greece, across different elevations and habitat types, we investigate the following research questions: 1. Is there a delay in the appearance of species and a progressive shortening of the flight period at higher elevations? We predict that phenological windows of activity will be shorter, better synchronized and delayed with elevation (Despland et al., 2012;Illán et al., 2012).
2. Are altitudinal patterns of butterfly phenology consistent with the temperature lapse rates recorded for each mountain system?
We expect to find steeper changes in emergence patterns when temperature lapse rate is steep and therefore climatic differences between elevations are more pronounced.
3. Do phenological patterns differ among different habitat types (agricultural areas, grassland, and forest)? We expect that phenology may vary with elevation at a different rate in different habitat types . For example, microclimates in forests that have a denser canopy compared to open habitats are less influenced by direct radiation (Scherrer & Körner, 2011), potentially leading to longer delays in emergence compared to open habitats (e.g., grasslands).

Do the responses of individual species follow consistent patterns
with elevation in terms of time of the appearance and duration of the flight period? We expect that univoltine species will show less pronounced altitudinal variation in phenology as a result of lesser adaptability compared with multivoltine species.

| Study system
Our study area consisted of two mountain regions that differ in geographic position, areal extent, biome, climate type, and topography. The Rodopi mountain chain (Rodopi hereafter: long. 24° 23′, lat. 41°23′; maximum elevation 2,323 m) is located in NE Greece, whereas Grammos (long. 20°50′, lat. 40°21′; maximum elevation 2,520 m) is located in NW Greece ( Figure S1, Table S1, but see also Zografou, Wilson, Halley, Tzirkalli and Kati (2017) for detailed descriptions). Both systems share a low human population density and associated low-intensity human activities, as well as high coverage by protected areas of the Natura 2000 network. The climate in Rodopi is at the transition between Mediterranean and a continental climate (Mavromatis, 1980) with a mean annual temperature of 11.4°C and mean annual precipitation of 1,200 mm, while the climate in Grammos is humid continental (Korakis, 2002) with a mean annual temperature of 8-12°C and mean annual rainfall of 1,500 mm.

| Butterfly sampling
Butterflies were recorded at 41 sites in Rodopi and 26 in Grammos.
The minimum distance between nearest neighboring sites was approximately 2 km (SD ± 0.5) so that each site effectively represents an independent sampling unit. The lack of spatial autocorrelation between nearby sites was verified in terms of alpha and beta components of diversity in a previous study where we investigated diversity patterns of butterflies and Orthoptera across different spatial scales (Zografou et al., 2017). Each mountain was partitioned into four elevation zones (0-500 m, 501-1,000 m, 1,001-1,500 m, and 1,501-2,000 m) and each zone contained sites representing the three dominant habitats found in the study system (agricultural fields, grasslands, and forests), with the exception that agricultural areas were not present above 1,500 m ( Figure S1, Table S1).
Permanent transect routes were established at two to six sites representative of each habitat type per altitudinal zone in each mountain range, recording geographic location (UTM) and elevation (m) using a hand-held GPS unit. On each site visit, the transect was walked at a steady pace under weather conditions that were suitable for butterfly activity (Pollard & Yates, 1993) recording all butterflies observed along a standardized length and width of 300 × 5 m.
Butterflies were captured with the help of hand net, identified in situ, and when necessary photographic material was also collected for confirming identification in the laboratory. We visited each site five times from April until August 2012 (Rodopi) and four times from May until August 2013 (Grammos-no sampling conducted in April due to unsuitable weather). Each transect was walked with a maximum sampling interval of 20 days between visits: This was the minimum interval which was feasible for a single field observer to achieve, given unpredictable weather and occasionally inaccessible sites particularly at higher elevations.

| Phenological descriptors
The timing and duration of flight periods were calculated to describe species' phenology along the altitudinal gradient. For each species, the timing of flight period was summarized per site as the weighted mean date (hereafter mean date) by summing counts per visit, ac- Mean date is a commonly used descriptor in phenological studies for butterflies and considered to be more reliable than other phenological measures such as the first day of adult appearance (Van Strien et al., 2008). In addition, as the occurrence of butterfly individuals in temperate species follows an approximately normal frequency distribution (Arce Crespo & Gutierrez, 2011), the use of mean date considers to be a safe approach (Moussus, Julliard, & Jiguet, 2010). We acknowledge the lack of multiple visits per month (e.g., weekly) but we emphasize that the main purpose is to examine relative differences in the degree of phenological shift, rather than to get unbiased estimates of the extent to which phenology change. We also calculated the duration of flight period as the standard deviation about the mean date (Brakefield, 1987). At the community level, we used all species for which the estimation of the two phenological descriptors could be generated (87) and when comparing the mean flight dates of butterfly communities between the two mountains and across different habitat types, only species present in both mountains (87)  the number of sites where the species was present, the minimum elevation, the maximum elevation, and the elevational range (maximum-minimum).

| Data analysis
For our first hypothesis, we investigated variation with elevation in the timing and duration of the flight period. We carried out linearmixed models where the mean date and standard deviation about the mean date were modeled as a function of altitude, mountain, and habitat. In addition, species were included as a random effect. Models were validated by checking for homoscedasticity and normality of the residuals (Zuur, Ieno, Walker, Saveliev, & Smith, 2009), and in all cases, diagnostic graphs showed that model assumptions were met ( Figure S2). For these models, altitude slope represented the delay (in days/km).
To investigate our second research question, we evaluated whether butterfly assemblages occurring in Rodopi have greater elevational delays in emergence compared with their counterparts in Grammos (considering species common to both mountain ranges), as a result of the different rates of climatic variation with elevation between the two mountains. We did this using standardized major axis (SMA) analysis. SMA is especially suitable when the prime interest is to inspect the slopes to see how each pair of variables is related to each other, rather than predicting Y (phenological descriptor) from X (elevation). In addition, SMA is a slope-fitting technique that shows how one variable scales against another, and slopes are fitted via a permutation test by minimizing the residual variance in X and Y dimensions simultaneously rather than Y alone (Domínguez et al., 2012;Falster & Westoby, 2005) resulting thus in a less biased outcome compared to traditional approaches such as ANCOVA (Warton, Wright, Falster, & Westoby, 2006). As a result, the sampling error which in our case is derived by the high topographic variability of mountain ranges can be minimized and biased slopes avoided (Legendre, 2001). Although SMA has been recommended particularly for allometric studies (Warton, Wright, & Wang, 2012), it can also be applied to ecological responses to environmental variables in the context of climate change (Zografou, 2015). Mountain was used as the grouping factor, and we discarded the first sampling in Rodopi (April 2012) in order to ensure that data (of four visits between May and August) were comparable between the two mountain ranges.
For the third research question, we used the same approach and tested whether butterfly assemblages that occur in forests have longer phenological delays (steeper slopes) with elevation compared to their counterparts in grasslands or agriculture areas.
To investigate our last research question regarding variation in the phenology of individual species with elevation, we ran general linear models for the 26 selected species, using the regression slope to estimate the delay in days/m. To account for between mountain and habitat variation, both terms were included in the models and p values were corrected using Benjamini and Hochberg (1995) adjustment method. We also tested whether the elevational delay was related to the number of sites where a species was present and the species' elevational range. Species that are present in more sites or species with wider elevational ranges may be expected to have longer delays compared to those whose distributions are limited to fewer sites or high elevations only, because the latter species may exhibit flight periods synchronized within a narrow phenological window, for example, avoiding the risk of unfavorable weather conditions in late summer (Illán et al., 2012).
To visualize our general linear-and mixed-effects models with partial residual plots, we extracted adjusted data using the "visreg" function in the "VISREG" package (Breheny & Burchett, 2017).
Finally, the "ggplot" function of the "ggplot2" package (Wickham, 2016) was used and ggplot2 library for the graphical representation of our results.

| Temperature lapse rate
We collected temperature data at each site using a Hobo data logger, to determine the gradient in seasonal temperature over elevation

| Phenological patterns at the community level
We found a positive and significant relationship between the variables "mean date" and "elevation" for the whole species pool

| Phenological patterns between mountains and across habitats
The relationship between mean date and elevation differed between mountains (LR test: 6.95, p = .007, n = 67) indicating a different rate at which butterfly assemblages delayed their appearance date with elevation ( Figure 2

| Phenological patterns at species level
We analyzed elevational patterns for 26 species: 20 had a positive slope when testing the relationship between the mean date and elevation, and six had a negative slope ( We found no interspecific evidence of effects on the species' elevational delay of the number of sites where species occurred (p = .98, n = 26) or the width of the species' elevational range (p = .40, n = 26).

| Temperature lapse rate
Considering both mountains and years, we found a significant de-

| Date of appearance
In the two Mediterranean mountains studied, flight dates occurred later for butterfly communities at higher elevations, in agreement with the few previous studies of Mediterranean mountain butterfly communities (de Arce Crespo & Gutiérrez, 2011;Gutiérrez & Menéndez, 1998;Illán et al., 2012). The flight dates of individual species also generally occurred later at higher elevations (see also Forister & Shapiro, 2003). On the basis of a temperature lapse rate of approximately 3°C per every kilometer in elevation increase, our findings suggest that a 1°C decrease in mean seasonal temperature could be associated with a 4.66-day phenological delay at the community level and an 7.71-day (average) phenological delay at the species level. A similar trend of a 3.7-day phenological delay for the entire butterfly community has been reported from Spain (de Arce Crespo & Gutiérrez, 2011;Illán et al., 2012).
The majority of univoltine species (70%) delayed the day of appearance with increasing elevation, whereas only 31.25% of the multivoltine species showed both delay and advance (Table 1).
Overall, shifts in phenology for less flexible species, as the ones with a single annual reproductive cycle, are less pronounced than showed here (Macgregor et al., 2019). It is not, however, the first time where univoltine butterflies seem to be more prone to develop adaptations against the dry and hot summer of Mediterranean region (Garcia- Barros, 1988). A potential butterfly strategy for increasing caterpillars' survival rate is to avoid the dry summer period and becoming more active on the cooler and wet months of early spring (Lopez-Villalta, 2010).
On the other hand, previous work suggests that species with multiple generations may take advantage of warming conditions by increasing the number of generations (Altermatt, 2010) and thus showing thermal plasticity in life cycle regulation (Van Dyck et al., 2015). Greater synchrony in time of emergence across temperature gradients for multivoltine species has also been interpreted as a possible sign of adaptation to local climatic conditions (Roy et al., 2015). The altitudinal delays we observe are thus likely to be subject to plastic variation depending on annual climatic conditions.
For example, Suggitt et al. (2012) found that butterflies occurring in  Illán et al., 2012), these rates are unlikely to represent fixed attributes of the species. seed production as the climate has warmed (Steltzer & Post, 2009).
Alternatively, negative species patterns might be regulated by an evolutionary adaptability to warmer climate, through an increased voltinism, despite the cooler local conditions at high elevation. An earlier appearance and prolonged flight period within areas above the timberline has also been reported for an alpine butterfly species

| Duration of the flight period
Almost a 3-day decline in the duration of the flight period for the community and an 8.86-day average decline considering responses by individual species over elevation are in agreement with previous findings in the Mediterranean area (de Arce Crespo & Gutiérrez, 2011;Illán et al., 2012). We argue, however, that the lack of significant individual responses for most of the species tested is not simply due to a lack of statistical power, given that interspecific variability in elevational delays did not appear to be associated with sample size or elevational range.
Difference in species' diet spectrum might drive phenological changes at the level of individual species. We noticed that out of the four species that shifted their flight periods along elevation, three were woody feeders and only one was herbaceous plant feeder.
According to Altermatt (2010), species feeding on herbaceous plants have smaller shifts in flight periods than the woody feeders because the second group has a narrow window to match the phenology to the flushing leaves. While herbs can produce leaves throughout the growing season, woody plants usually flush their leaves simultaneously (Feeny, 1976), and only for a short time of the year can accommodate the needs of herbivores for fresh and palatable resources.
For the four species that exhibited a negative pattern in their duration of the flight period, three were multivoltine species and one was a univoltine species. In this context, the more pronounced phenological patterns for multivoltine life cycles are indeed a function of the species' voltinism.

| Phenological patterns between mountains
The most striking feature is the inconsistency between elevational delays of butterfly assemblages between the two mountains with respect to temperature lapse rate.  (Korakis, 2002;Xirouchakis, 2005). suggesting no such effect.
We argue that the steeper temperature lapse rate (5°C) along the elevation in Grammos may have driven species to better synchronize their activity resulting in a smaller delay overall. Empirical evidence suggests that populations from more variable environments have higher levels of plasticity which could preadapt them to extremes (Chevin & Hoffmann, 2017). When such extremes are lacking, it is logical to assume that species responses are not masked by phenotypic plasticity and therefore are steeper and more pronounced. Another explanation could be that species in Rodopi are closer to their upper thermal limits, and it is unlikely to evolve physiological tolerances to increased temperature (Araújo Miguel et al., 2013;Mills Simon et al., 2017). Note: Int: intercept, SE: standard error, subscript md corresponds to mean date and d to the duration of the flight period. Univoltine species are indicated by the subscript letter u, while the rest have more than one generation. Significance codes: 0 "***" 0.001 "**" 0.01 "*" 0.05 "ns" nonsignificant.
The species are in alphabetical order. The number of sites occupied, the minimum and maximum elevations (m), and the elevational range for each species are included too.
a Selected species are species recorded in more than three sites with at least two records per site and species overwintering as egg, pupae, or larvae, excluding thus early spring flyers and species overwintering as adults for which phenology may not be recorded comprehensively.
permanent meteorological stations within each region might help to inform our findings.

| Predictions and conclusions
On the basis of the different climatic scenarios proposed for the Eastern Mediterranean and Middle East, the mean temperature rise will be about 1-3°C in the near future (2010-2039), 3-5°C by midcentury (2040-2069), and 3.5-7°C by the end of the century (2070-2099) (Lelieveld et al., 2013). given the complexity and dynamism of the natural system, such radical changes of temperature are likely also to change many other contributing factors too, such as weather conditions at different times of year, as well as ecological community structure, where we are likely to see warm-adapted species expanding at the expense of cold-adapted ones (Zografou et al., 2014). A further aspect of system complexity is the ongoing forest encroachment that tends to counteract climate change, benefiting woodland species at the expense of others (Slancarova et al., 2016).
While it is difficult to foresee how organisms are going to cope under the ongoing changes in climate, it is possible that advanced emergences could threaten serious trophic disruption between interacting groups. Our findings both confirm an earlier and prolonged activity at lower elevations overall. At the same time, confound expectations for ectotherms such as signs of earlier appearance in high elevations for multivoltine organisms and more pronounced shifts in flight periods for the woody feeders challenge the idea that these species assemblages have special thermal traits that confer adaptive advantage under new conditions. National Park for research support.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
VK, JMH, RJW, AG, and KZ conceived and designed the experiments. KZ performed the experiments, and KZ and GCA analyzed the data. KZ wrote the first draft of the manuscript, and AG, RJW, JMH, VK, and GCA contributed to subsequent versions of the article and agreed on the final version to be published.

O PE N R E S E A RCH BA D G E
This article has earned an Open Data Badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at https ://datad ryad.org/stash/ share/ 0wyAa Q Clpd viKiC 3L8Ufg_U0pOc HfHVX 8PwK9 tUfKI0.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data connected to the manuscript with the title "Butterfly phenology in Mediterranean mountains using space-for-time substitution" have been deposited in a publicly accessible repository Dryad: https ://datad ryad.org/stash/ share/ 0wyAa QClpd viKiC 3L8Ufg_ U0pOc HfHVX 8PwK9 tUfKI0.