Phenological sensitivity and seasonal royalsocietypublishing.org/journal/rspb variability explain climate-driven trends in Mediterranean butterflies Pau Colom1, Miquel Ninyerola2, Xavier Pons3, Anna Traveset1 and Research Constantí Stefanescu4,5 Cite this article: Colom P, Ninyerola M, Pons 1 Global Change Research Group, Institut Mediterrani d’Estudis Avançats (IMEDEA-CSIC-UIB), Miquel Marqués 21, X, Traveset A, Stefanescu C. 2022 Phenological 07190 Esporles, Mallorca, Balearic Islands, Spain 2 sensitivity and seasonal variability explain Grumets Research Group, Departament de Biologia Animal, Biologia Vegetal i Ecologia. Edifici C. Universitat Autònoma de Barcelona, 08193 (Bellaterra, Barcelona), Catalonia, Spain climate-driven trends in Mediterranean 3 Grumets Research Group, Departament de Geografia. Edifici B, Universitat Autònoma de Barcelona, butterflies. Proc. R. Soc. B 289: 20220251. 08193 (Bellaterra, Barcelona), Catalonia, Spain 4 https://doi.org/10.1098/rspb.2022.0251 Natural Sciences Museum of Granollers, Francesc Macià 51, 08402 (Granollers, Barcelona), Catalonia, Spain 5 Centre de Recerca Ecològica i Aplicacions Forestals (CREAF-CSIC-UAB), Universitat Autònoma de Barcelona, 08193 (Cerdanyola de Vallès, Barcelona), Catalonia, Spain PC, 0000-0003-0309-8886; CS, 0000-0001-8952-7869 Received: 8 February 2022 Accepted: 28 March 2022 Although climate-driven phenological shifts have been documented for many taxa across the globe, we still lack knowledge of the consequences they have on populations. Here, we used a comprehensive database com- prising 553 populations of 51 species of north-western Mediterranean butterflies to investigate the relationship between phenology and population Subject Category: trends in a 26-year period. Phenological trends and sensitivity to climate, Ecology along with various species traits, were used to predict abundance trends. Key ecological traits accounted for a general decline of more than half of Subject Areas: the species, most of which, surprisingly, did not change their phenology ecology under a climate warming scenario. However, this was related to the regional cooling in a short temporal window that includes late winter and early spring, during which most species concentrate their development. Finally, Keywords: we demonstrate that phenological sensitivity—but not phenological global warming, insect phenology, trends—predicted population trends, and argue that species that best population trends, species’ life-history traits, adjust their phenology to inter-annual climate variability are more likely host plant specialization, long-term trends to maintain a synchronization with trophic resources, thereby mitigating possible negative effects of climate change. Our results reflect the importance of assessing not only species’ trends over time but also species’ abilities to respond to a changing climate based on their sensitivity to temperature. Authors for correspondence: Pau Colom e-mail: pcolom@imedea.uib-csic.es Constantí Stefanescu e-mail: cstefanescu@mcng.cat 1. Background The effects of anthropogenic climate change have been documented for a wide range of taxa [1], of which insects—and, specifically, butterflies—are one of the best-studied groups ever since research into this question began [2,3]. The remarkable sensitivity of these ectothermic organisms to changes in biotic and abiotic conditions, together with the availability of long-term databases covering large geographical areas generated by citizen science projects, have made this group a study model [4]. The ecological literature in recent decades has highlighted how changes in the population size [5–7], distribution [3,8,9] and phenology [5,10,11] of butterflies are being driven by an ever-warmer climate. Nevertheless, how population trends are being affected by these various changes is less clear (but see [12,13]). The advancement of spring emergence due to increasing temperatures is the Electronic supplementary material is available most frequently recorded phenological trend in butterflies [5,14,15]. However, online at https://doi.org/10.6084/m9.figshare. c.5954002. © 2022 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited. intra- and interspecific variation in the degree of this response to avoid unreliable flight curves. A generalized additive model 2 may vary due to latitudinal [16,17] and seasonal [15,18] (GAM) was then fitted to the site counts to generate the flight curve of a species for the site and year (see below, ‘generation royalsocietypublishing.org/journal/rspb climatic variation. In the Mediterranean basin, in particular, there have been strong seasonal differences in the magnitude of variables’), but excluded cases in which the model did not converge or failed to find at least one peak of abundance. Finally, of warming [19]. Therefore, we would expect the strength of we discarded (i) migratory species (i.e. when records mostly cor- phenological trends to depend on the specific period of the respond to immigrant individuals; e.g. Vanessa cardui); (ii) multi- year in which the development of the immature stages is voltine species in which the different generations greatly overlap affected by the ambient temperature. However, to date, no (e.g. Pararge aegeria and Polyommatus icarus) and (iii) species for study has yet assessed this hypothesis. which we had data from less than three populations. Our final Despite the considerable amount of data documenting phe- dataset included 553 populations of 51 species that were statisti- nological shifts, little is known about how changes in the cally robust enough to estimate long-term trends in both timing of emergence influence butterfly population dynamics. abundances and phenologies. Advancing emergence may have several positive implications. First, it can reduce the time the immature stages need to (b) Ecological and life-history species traits Proc. R. Soc. B 289: 20220251 develop, and thus, the time they are exposed to potential pre- We used six traits that could potentially explain trends in abun- dators and parasitoids [20,21]. Second, it may increase the dance and phenology, and differences in species’ sensitivity to time available to complete generations, with potential benefits temperature [11,12,29,30]. Three traits were categorical: voltinism for abundance trends [22,23]. Finally, it may ensure good syn- (univoltine versus multi-voltine), overwintering stage (egg, chronization with the phenology of plants used as nectar larvae, pupa or adult) and larval diet composition (grass, forb or sources, which are also expected to advance their flowering tree); the other three traits were continuous variables related to due to warming [24]. On the other hand, some studies have adult habitat preference and specialization, and to larval trophic suggested that phenological advances could in fact be mala- specialization: (i) a habitat species specialization index (SSI) [31]; daptive. For instance, they could be detrimental if they lead (ii) an index to estimate butterfly preferences for open or closed to the addition of new generations that do not have enough habitats (TAO index, [32]) and (iii) the host plant index (HPI, time to enter diapause before the arrival of winter (i.e. develop- [33]) that quantifies the trophic specialization of a butterfly species in the larval stage based on the number of families, genera and mental trap; [23,25]). Taken together, all these possible positive plant species used in the study region. This last variable was nor- and negative effects involve great complexity to produce accu- malized using a square-root transformation. See electronic rate predictions regarding the effects of changing the time of supplementary material, appendix, text S1 for further details of emergence on abundance trends. the trait calculations. Recent research has shown positive effects of phenological advancement on abundance trends in Lepidoptera of North and Central Europe [12,13]. However, phenological effects on (c) Temperature data abundance have only been investigated using phenological Month-by-month mean air temperature (minimum, mean and trend as the predictor variable, while none of these studies maximum) raster surfaces were obtained using a two-step used phenological sensitivity to climate as a predictor of popu- (global and local) spatial interpolation (see [34] for details). First, a multiple linear regression using geographical predictors (alti- lation dynamics. Here, we explored the idea that species that tude, latitude, continentality and solar radiation) was performed show high plasticity of phenology to inter-annual climate for each month and climatic element. Secondly, an inverse distance variability have a greater ability to cope with climate change, weighting residual local interpolation was used to correct these as this trait will favour a greater synchronization with trophic surfaces. This method is based on data from ca 3500 meteorological resources that can be highly plastic (e.g. [26]), and therefore, stations run by the State Meteorology Agency (AEMET) and a it would increase individual fitness. We hypothesize that quasi-automatic quality control suite, applied to deal with this species with high sensitivity to climate—that is, species more massive dataset and to remove the effects of non-climatic factors likely to advance their phenology in warmer years but also (i.e. short length series, non-stable or poorly curated data, or to delay it in colder years—have experienced more positive station relocation accompanied with geolocation refinements). or fewer negative abundance trends in recent decades. The obtained maps gave a cross-validated root mean-squared error of ca 1°C for all months. Initial 100 m surfaces were densified into 2 m surfaces (using bilinear interpolation) to create new ver- tices in the butterfly sampling sites. Then, we enriched these 2. Material and methods vertices with the pixel values from the climate rasters and com- puted the mean temperatures for each site. (a) Study area and butterfly data The study area lies in the north-west Mediterranean basin (ca 32 600 km2) in an area in which butterflies have been systematically (d) Abundance and phenological estimates monitored since 1994 within the framework of the Catalan and For each population, we calculated annual abundance and Andorran butterfly monitoring schemes (CBMS and BMSAnd; phenological estimates. Annual abundance was calculated as the see www.catalanbms.org) (electronic supplementary material, logarithm of the total number of individuals recorded in a year. appendix, figure S1). We divided the total number of individuals by the number of record- Butterfly abundance data were recorded weekly from March ing events as in some years, the transects were not sampled during all to September (i.e. ca 30 recording events each year depending on the 30 weeks due to bad weather conditions in a few weeks. To esti- the weather conditions) along fixed transects following a stan- mate population phenology, we followed the approach used by dard methodology [27,28]. We used data from a 26-year period Macgregor et al. [12] and fitted abundance data using GAM (1994–2019) from 59 sampling sites that had been monitored models with a Poisson error distribution with a restricted maximum for at least 10 years. For all possible combinations of species, likelihood. This method is less sensitive to variation in sample sizes sites and years, we restricted our dataset to cases in which the and species detectability and provides more accurate estimates than focal species were recorded during at least five weeks at a site methods based on first appearance dates [35]. We used the date (Julian day) of the first annual emergence peak predicted by the (a) 3 models as the annual phenological estimate for both univoltine 20 annual mean temperature (°C) and multi-voltine species (see electronic supplementary material, royalsocietypublishing.org/journal/rspb appendix, text S2 for more details and specifications on how these 18 calculations were refined for each species). 16 14 (e) Data analysis 12 (i) Temperature trends We first analysed the overall temperature trends for our 59-site 10 dataset. We carried out generalized linear mixed models 8 (GLMM) for the mean monthly and annual temperatures 1994 1999 2004 2009 2014 2019 (response variables), using year as a predictor and site as a random effect, and assuming a Gaussian error distribution. We (b) 2 repeated the same analysis with maximum and minimum temp- trend of mean temperature Proc. R. Soc. B 289: 20220251 eratures as response variables. We also calculated the mean annual trends for each site separately with linear regressions. Finally, we tested the spatial variability of our dataset using a cat- 1 egorization of the sites’ climatic region with a three-level factor (alpine and subalpine, Mediterranean mesic and Mediterranean xeric). We performed a generalized linear model (GLM) with 0 the slope of the site trend (i.e. annual mean temperature versus year) as the response variable and the climate region as the pre- dictor variable, again adjusting errors to a Gaussian distribution. –1 Jan Feb Mar Apr May June July Aug Sep Oct Nov Dec (ii) Abundance and phenological trends Figure 1. Annual and monthly mean temperature trends over a 26-year Abundance and phenological trends were calculated at both popu- period (1994–2019). (a) Positive significant trend in annual mean tempera- lation (n = 553) and species (n = 51) levels. Population trends were ture. The red line shows the linear regression line over time. Bars represent calculated as the slope of a linear regression with abundance (i.e. the standard deviation between the 59 study sites. (b) Violin plots show the logarithm of the mean abundance) and phenology (i.e. Julian frequency distribution of the slopes of the linear models between tempera- day) as response variables and year as the predictor variable, with errors adjusted to a Gaussian distribution. Species trends ture and year for the 59 study sites in each month. The black dashed line in were calculated with the same structure but using GLMMs and each violin plots represents the median of the distribution. The red dotted controlling for the site factor as a random effect. line at the zero value on the y-axis is shown to distinguish between negative and positive values. (Online version in colour.) (iii) Phenological sensitivity To assess the sensitivity of species phenology to ambient temp- this, we excluded three species as outliers because they have posi- erature, we used a time-window approach [15]. We calculated tive values of phenological sensitivity (i.e. responded with the mean temperatures of 36 potential time-windows ranging phenological delays to an increase in temperature; see electronic from one to three months from September of the previous year supplementary material, appendix, text S3) in contrast with the to September of the current year. We used these potential time- other 48 species with negative values (i.e. responded with expected windows to identify the period of the year with the highest temp- phenological advancements to an increase in temperature) erature effect on the phenology of each species (hereafter, the (figure 2c). We used the absolute value of the slope of the trend critical period). For each time-window and species, we fitted the over time (hereafter, the absolute phenological trend), which quan- phenological data using a GLMM with temperature as a predic- tifies the degree of change without considering its direction (i.e. tor, site as a random effect and a Gaussian error distribution. For either advancement or delay). This variable was first normalized each species, we selected the best-fit model according to the with a square-root transformation because of its skewed distri- Akaike information criteria. We used the slope of the best-fit bution. Finally, we tested the effect of phenological sensitivity, models to estimate species’ phenological sensitivity. phenological trend and absolute phenological trend on abundance trends. We used both the real and absolute values of the trends as variables to test separately the effect of the direction and the mag- (iv) Species trends related to ecological and life-history traits nitude, respectively, of the changes in the emergence peak on the To test the effects of species traits on phenological sensitivity abundance trends. We first ran the models including only a and trends in abundance and phenology, we used phylogenetic single predictor variable (either phenological sensitivity, phenolo- generalized least-squares (PGLS) to control for the potential phylo- gical trend or absolute phenological trend) and then also included genetic signal, as we assumed that more closely related species are the significant species traits related to abundance trends. more likely to have similar ecological and life-history traits [36]. To All analyses were carried out using R v. 4.1.1 [38]. The fol- construct the phylogenetic tree of our species subset, we used lowing packages were used: mgcv [39] to fit GAM models, lme4 recently published data on the phylogeny of European butterflies [40] to conduct GLMMs, ape [41] to carry out the PGLS analyses, [37] and PGLS to analyse the effect of six traits (see above, ‘Ecologi- and MuMIn [42] for model comparison. cal and life-history species traits’) on (i) phenological sensitivity, (ii) phenological trends and (iii) abundance trends. 3. Results (v) Relationship between abundance trends and phenological variables (a) Trends in local temperature We performed several PGLS models to test to what extent phenol- Our models estimated an increase of 0.59°C in mean tempera- ogy can explain abundance trends. First, we tested whether greater tures (t = 14.02, p < 0.01) (figure 1), 0.55°C in maximum phenological sensitivity results in a steeper phenological trend. For temperatures (t = 8.64, p < 0.01) and 0.42°C in minimum (a) but decreased significantly in February, March and May 4 15 e (figure 1; electronic supplementary material, appendix, table e as cr royalsocietypublishing.org/journal/rspb in S2). A similar seasonal pattern was observed when we used decline stable an independent dataset of average temperatures recorded in 10 a network of 24 automatic weather stations for the entire no. species study region (Meteorological Service of Catalonia, data not shown). No differences in temperature trends were observed 5 between climatic regions, which allowed us to assume that spatial variability is not significant in our dataset. 0 (b) Trends in abundance and phenology –0.030 –0.020 –0.010 0 0.010 0.020 0.030 Of the 553 analysed populations, 33% showed significant abundance trend trends in abundance, of which 75% were negative (i.e. popu- (b) lations declined during the study period). Similar results Proc. R. Soc. B 289: 20220251 . 20 were observed at species level, with a tendency for abundances anc adv to decline (figure 2a); 53% of species significantly declined in delay 15 abundance, while only 12% showed significant positive stable trends (electronic supplementary material, appendix, table no. species S3). Regarding phenology, only 3% of the populations 10 advanced their emergence peaks, while 3% delayed them. At species level, 8% and 15% of species advanced and delayed 5 their peaks, respectively (figure 2b; electronic supplementary material, appendix, table S3). 0 –0.6 –0.4 –0.2 0 0.2 0.4 0.6 0.8 1.0 1.2 (c) Butterfly sensitivity to climate phenology trend The effect of temperature on butterfly emergence peak was (c) negative for most species, indicating a clear trend to advance delay 15 the flight period as temperatures increase (figure 2c). Critical ble sta periods varied between species (figure 3). Temperature effects on emergence peaks were greatest in the autumn of 10 advancement the preceding year (for Callophrys rubi and Vanessa atalanta), no. species in winter (for most species) or in spring–summer (e.g. for Favonius quercus and Hesperia comma). Up to 86% of the 5 species included January in their critical period, 57% Febru- ary and 22% March. An increase in the annual temperature of the critical period produced significant advances in 90% 0 of species’ phenology (electronic supplementary material, –0.7 –0.6 –0.5 –0.4 –0.4 –0.2 –0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 appendix, table S3). Only two species responded to warming phenological sensitivity with a significant phenological delay (Favonius quercus and Hipparchia fagi; see electronic supplementary material, appen- Figure 2. Abundance, phenology and sensitivity estimates for the 51 studied dix, text S3), while three species did not respond significantly butterfly species. The estimates are the slopes of the GLMMs with (a) the to temperature changes (Libythea celtis, Hipparchia semele and logarithm of mean abundance as response variable and year as predictor; Spialia sertorius). (b) Julian day of the emergence peak as response variable and year as pre- dictor; (c) Julian day of the emergence peak as response variable and mean temperature of the critical period as predictor. The dashed line at the zero (d) Interspecific variability in butterfly trends value on the x-axis is shown to distinguish between negative and positive values. Negative values of phenological trends and phenological sensitivity modulated by species traits Several traits showed to be significant predictors of abun- indicate advances in butterfly emergence in time or responses to increasing dance trends, including HPI, voltinism, SSI and TAO index temperature, respectively, while positive values indicate delays in butterfly (figure 4; electronic supplementary material, appendix, emergence. Pie charts show the percentage of species in each category. table S4). Species using fewer host plants (i.e. trophic special- (Online version in colour.) ists), univoltine species, habitat specialists and species with a preference for more open habitats had the greatest declining temperatures (t = 7.18, p < 0.01), on average, for the 59 sites trends. However, the model comparison indicated that during the 26-year period (1994–2019) (electronic supplemen- larval trophic specialization (measured as HPI) was the best tary material, appendix, figure S2). According to individual predictor explaining abundance trends. Larval trophic linear models, 51% of sites showed a significant increase in specialization was also the only significant predictor of phe- mean temperature but no site showed a significant decrease nological trends, which suggests that specialized species in this variable (electronic supplementary material, appendix, experience a greater degree of phenological delay (electronic table S1). Nevertheless, monthly temperature trends for supplementary material, appendix, table S4). Species sensi- all sites showed remarkable seasonal variability. Mean temp- tivity was not significantly predicted by any species trait eratures increased significantly in nine months of the year (electronic supplementary material, appendix, table S4). 5 Callophrys rubi Vanessa atalanta royalsocietypublishing.org/journal/rspb Coenonympha dorus Plebejus argus Lysandra bellargus Lysandra coridon Coenonympha arcania Gonepteryx cleopatra Pyronia cecilia Zerynthia rumina Anthocharis cardamines Pseudophilotes panoptes Anthocharis euphenoides Colias crocea Euchloe crameri Leptidea sinapis Pieris brassicae Proc. R. Soc. B 289: 20220251 Pieris napi Coenonympha pamphilus Aporia crataegi Aphantopus hyperantus Melitaea didyma Melitaea cinxia Eurphydryas aurinia Aglais urticae Speyeria aglaja Iphiclides feisthamelii Thymelicus sylvestris Thymelicus acteon Spialia sertorius Pyrgus malvoides Ochlodes sylvanus Pyronia bathseba Melanargia lachesis Maniola jurtina Lasiommata megera Hipparchia hermione Coenonympha arcania Hipparchia fidia Libythea celtis* Satyrium esculi Argynnis paphia Charaxes jasius* Gonepteryx rhamni Brintesia circe Hipparchia semele Hipparchia statilinus Celastrina argiolus* Favonius quercus Hesperia comma Hipparchia fagi Sep Oct Nov Dec Jan Feb Mar Apr May June July Aug Figure 3. Species’ critical periods. The period from September of the previous year to September of the current year in which the sensitivity (i.e. the relationship between the temperature and the emergence peak) is greatest. Lines represent the critical period of each species: in red, months with positive trends in temperature; in blue, months with negative trends in temperature (figure 1; electronic supplementary material, appendix, table S2). Red and blue asterisks show respectively positive and negative significant trends of the critical period. Points represent the mean peak day of each species. (*) We used the peak inactivation of hibernating adults for Libythea celtis and the peak in the second generation for Charaxes jasius and Celastrina argiolus (electronic supplementary material, appendix, text S2). (Online version in colour.) (e) Effects of phenological responses to climate on more positive or fewer negative abundance trends (figure 4). Neither the absolute phenological trend nor the real phenologi- abundance trends cal trend variables had any significant effect on abundance The absolute phenological trend was significantly associated trends, regardless of including the significant species traits in with phenological sensitivity (electronic supplementary the model. material, appendix, table S5). The greater the sensitivity to temperature a species had, the greater change in the emergence timing experienced ( phenology delay or advancement) during the study period. Phenological sensitivity also significantly 4. Discussion predicted abundance trends, both as the sole predictor and when all significant species traits (electronic supplementary (a) Phenological trends material, appendix, table S5) were included in the model. Overall, our findings contrast with the expected pattern of Species with greater sensitivity to temperature experienced phenological advancement under a global warming scenario 6 phenology ecological traits life-history traits royalsocietypublishing.org/journal/rspb 0.02 0.02 0.06 p = 0.688 p < 0.001 p = 0.057 0.04 0.01 0.01 abundance trend abundance trend abundance trend 0.02 0 0 0 –0.02 –0.01 –0.01 –0.04 n = 4 n = 30 n = 12 n = 5 –0.02 –0.02 –0.06 –0.6 –0.4 –0.2 0 0.2 0.4 0.6 0.2 0.4 0.6 0.8 1.0 egg larva pupa adult phenological trend host-plant index (HPI) overwintering stage 0.02 0.02 0.06 p = 0.967 p < 0.001 p = 0.513 0.04 0.01 0.01 abundance trend abundance trend abundance trend 0.02 Proc. R. Soc. B 289: 20220251 0 0 0 –0.02 –0.01 –0.01 –0.04 n = 20 n = 23 n = 8 –0.02 –0.02 –0.06 0 0.2 0.4 0.6 0.8 0.6 0.8 1.0 1.2 1.4 1.6 grass forb tree absolute phenological trend species specialization index (SSI) larval diet composition 0.02 0.02 0.06 p = 0.048 p = 0.032 p = 0.041 0.04 0.01 0.01 abundance trend abundance trend abundance trend 0.02 0 0 0 –0.02 –0.01 –0.01 –0.04 n = 16 n = 35 –0.02 –0.02 –0.06 0.2 0.3 0.4 0.5 0.6 0.7 0.2 0.4 0.6 0.8 univoltine multivoltine phenological sensitivity TAO index voltinism Figure 4. Predictor effect plots for abundance trends from PGLS tests. Black points and violin plots show the distribution of the data. Blue lines and points depict model-predicted relationships between abundance and independent variables. Dashed blue lines and bars indicate 95% confidence intervals. The dashed red line at the zero value on the y-axis is shown to distinguish between negative and positive values. (Online version in colour.) [1,5]. Indeed, recent work based on similar data and lengths pattern, as these specialists depend on very different types of time series to our study supports the idea that in recent of plants (e.g. grasses, forbs and trees) and include species decades, phenological advancements may not have been as with a wide range of emerging times over the season. More common as predicted [13]. Only 8% of the species in our work is needed to understand this relationship, which dataset underwent significant advances and, surprisingly, could have conservation implications if increasing delays 15% of the species in fact delayed their phenology. Reports eventually become maladaptive for these species. of phenological delays often contain sampling errors (i.e. short time series; [43]). However, this does not seem to be the case here because we used time series with an average (b) Abundance trends length of 24.9 years for our calculations. Although mean Our analyses once again reveal a worrying decline in butter- temperatures in the region increased by 0.59°C, the decrease fly populations in the Mediterranean basin [29], which is in mean temperatures in a few particular months could be even more severe than in other regions of Europe [45]. key to explain the variability in species’ phenological More than half of our study species decreased in abundance, trends. In fact, our results are consistent with the seasonal and probably some of the rarest and most specialized species variability of climate trends described for the Mediterranean excluded from our analyses are those undergoing the most region, with a greater contribution to warming in summer serious declines. We identified several species traits related and autumn than in winter and early spring [19,44]. The to this pattern. Habitat specialists declined more than gener- observed negative trends of mean temperatures in February, alists as they are more affected by changes in land use (e.g. March and May could offset the general warming trend [46]). Specifically, species preferring open habitats had greater and also explain the delays in some species, as these negative trends, probably due to land abandonment and veg- months are part of the critical period of development in the etation encroachment in habitats [32]. Contrary to results of a immature stages of most species. As expected, the most sen- previous study [29], we found that univoltine species under- sitive species to temperature were those experiencing the went steeper declines than multi-voltine species, which adds steeper phenological trends in absolute terms, thereby con- to the current debate in the literature: the greatest declines firming the hypothesis that butterfly phenology is driven have been reported for univoltine species [12,47,48] but also mainly by climate. Interestingly, species with greater trophic for multi-voltine species [29,49]. Interestingly, the best predic- specialization in their larval stages delayed their phenology tor of abundance trends was larval trophic specialization (i.e. more than generalists. We have no explanation for this HPI index). Larval specialists, like habitat specialists, are usually the most affected by habitat loss and degradation. In therefore minimize the problems arising from the loss of 7 addition, larval specialists are the species that are most poten- synchrony with these resources (see, e.g. [24]). However, royalsocietypublishing.org/journal/rspb tially affected by a loss of phenological synchrony with their although these explanations seem plausible, further empirical host plants due to climate change [50]. Indeed, larval special- work is needed to investigate whether phenological plasticity ists were the group of species delaying most, or advancing actually increases butterfly synchrony with larval host plants the least, their phenology, which could have led to a greater and adult nectar plants. Further research will also be needed mismatch with their host plants. to identify the traits explaining phenological sensitivity and thus facilitate the identification of vulnerable species to climate change. (c) Effects of phenological shifts on abundance Theoretically, climate-driven changes in phenology can have Abundance trends were found to be associated with phenolo- both positive and negative demographic consequences. Shifting gical sensitivity, but not to phenological trends, confirming phenology in response to climate can be maladaptive under the hypothesis that at least for Mediterranean butterflies, phe- some circumstances (e.g. [25]), but on the other hand, in many nological sensitivity is a better predictor than phenological cases, it will improve the performance of individuals by increas- trends per se for abundance trends. The effect of phenological Proc. R. Soc. B 289: 20220251 ing the synchronization with other trophic levels. Here, we show sensitivity on abundance trends is robust since it remained that greater phenological sensitivity can prevent population significant when we controlled for the phylogenetic signal decline in butterflies and, probably, in other ectotherm organ- and for the effect of the ecological and life-history traits isms. Given that phenological climate-driven trends have intrinsically related to abundance trends. In the UK, Macgre- been documented in many taxa across a high variety of ecosys- gor et al. [12] found a positive effect for phenological tems and habitats, it is essential to better understand their advancement on abundance trends, but as in our work, this demographic consequences. We call for further research to effect has not been found in other regions [13,30]. The effect disentangle the ecological consequences of phenological in the UK was only found for multi-voltine species, so the climate-driven trends, assessing not only species’ trends over relatively small number of multi-voltine species in our time but also species’ sensitivity to temperature and their ability study dataset after removing the multi-voltine species with to respond to a changing climate. overlapping generations (n = 16) in comparison to the dataset of Macgregor et al. [12] (n = 39) could have prevented us finding the same pattern. However, multi-voltine species Data accessibility. Data are available for the reviewers and editors in the electronic supplementary material [55] and at the following link: have a greater sensitivity to increase their flight duration https://doi.org/10.6084/m9.figshare.19137512. than univoltine species [30], and therefore, this pattern Authors’ contributions. P.C.: conceptualization, formal analysis, investi- could be more related to the positive effects of lengthening gation, methodology and writing—original draft; M.N.: data the flight period [49] and the increases in voltinism described curation, methodology and writing—original draft; X.P.: data cura- in northern latitudes (e.g. [22]). Other strategies not analysed tion, funding acquisition and methodology; A.T.: funding in our study, including a northward shift in range bound- acquisition, supervision and writing—review and editing; C.S.: con- ceptualization, supervision, writing—original draft and writing— aries, could also be associated with positive abundance review and editing. trends [13], although this particular strategy will always be All authors gave final approval for publication and agreed to be highly dependent on suitable habitat availability [51], and it held accountable for the work performed therein. is possible that only species that already have stable or posi- Conflict of interest declaration. We declare we have no competing interests. tive trends could actually successfully undertake range Funding. The CBMS is funded by the Departament d’Acció Climàtica, expansions [52]. Alimentació i Agenda Rural de la Generalitat de Catalunya and the Although it is often argued that phenological sensitivity Diputació de Barcelona, and the BMSAnd by the Andorran Govern- ment. P.C. is funded by a PhD fellowship financed by the Govern de to changing climatic conditions may increase or maintain les Illes Balears (grant no. FPI-CAIB-2018). This study was partially population abundance, this has rarely been tested [53]. Our funded by the Catalan Government under grant no. SGR2017-1690 results actually suggest that phenological sensitivity could to X.P., by the Spanish MCIU Ministry through the NEWFORLAND be an adaptation to the inter-annual variation of weather research project (grant no. RTI2018-099397-B-C21 MCIU/AEI/ERDF, conditions in order to maintain synchrony with resources, EU) to X.P. and the DEPICT research project (grant no. PID2020- 114324GB-C2) to A.T. which have been shown to respond to the increase of temp- eratures by advancing their phenology ([26,54] and Acknowledgements. We are very thankful to all the CBMS and BMSAnd volunteers for collecting the butterfly data. We would also like to references therein). This could explain why less sensitive thank Andreu Ubach and Ferran Páramo for helping with the data species had sharper declines than highly sensitive species, management, Víctor Lecegui for calculating the Host Plant Index, and why they could potentially be the most negatively Antoni Palerm and Joshua Borràs for assisting in the construction affected in the current context of climate change. Further- of the loops in R for data analysis, Clàudia Pla-Narbona for preparing electronic supplementary material, appendix, figure S1, Marc more, the positive relationship between the degree of larval Prohom for his useful comments on climate trends in Catalonia, trophic specialization and population trends could also give and Michael Lockwood for revising the English version. Climate indirect support to this hypothesis, as it is the more generalist data were provided by the State Meteorological Agency (Ministry species that can find resources with different phenologies and for the Ecological Transition and the Demographic Challenge). References 1. Parmesan C, Yohe G. 2003 A globally coherent systems. Nature 421, 37–42. (doi:10.1038/ 2. Dennis RL. 1993 Butterflies and climate change. fingerprint of climate change impacts across natural nature01286) Manchester, UK: Manchester University Press. 3. Parmesan C et al. 1999 Poleward shifts in 13, 1860–1872. (doi:10.1111/j.1365-2486.2007. butterfly community. Sci. Rep. 11, 1–14. (doi:10. 8 geographical ranges of butterfly species associated 01404.x) 1038/s41598-021-82723-1) royalsocietypublishing.org/journal/rspb with regional warming. Nature 399, 579–583. 17. Post E, Steinman BA, Mann ME. 2018 Acceleration 31. Stefanescu C, Carnicer J, Peñuelas J. 2011 (doi:10.1038/21181) of phenological advance and warming with latitude Determinants of species richness in generalist and 4. Parmesan C. 2003 Butterflies as bioindicators for over the past century. Sci. Rep. 8, 1–8. (doi:10. specialist Mediterranean butterflies: the negative climate change effects. In Butterflies: ecology and 1038/s41598-018-22258-0) synergistic forces of climate and habitat change. evolution taking flight (eds CL Boggs, WB Watt, 18. Colom P, Traveset A, Carreras D, Stefanescu C. 2021 Ecography 34, 353–363. (doi:10.1111/j.1600-0587. PR Ehrlich). Chicago, IL: University of Chicago Press. Spatio-temporal responses of butterflies to global 2010.06264.x) (doi:10.7208/9780226063195) warming on a Mediterranean island over two 32. Ubach A, Páramo F, Gutiérrez C, Stefanescu C. 2020 5. Roy DB, Sparks TH. 2000 Phenology of British decades. Ecol. Entomol. 46, 262–272. (doi:10.1111/ Vegetation encroachment drives changes in the butterflies and climate change. Glob. Change Biol. 6, een.12958) composition of butterfly assemblages and species 407–416. (doi:10.1046/j.1365-2486.2000.00322.x) 19. Bartolini G, Di Stefano V, Maracchi G, Orlandini S. loss in Mediterranean ecosystems. Insect Conserv. 6. Forister ML, McCall AC, Sanders NJ, Fordyce JA, 2012 Mediterranean warming is especially due to Diver. 13, 151–161. (doi:10.1111/icad.12397) Thorne JH, O’Brien J, Waetjen DP, Shapiro AM. 2010 summer season. Theor. Appl. Climatol. 107, 33. Middleton-Welling J et al. 2020 A new Proc. R. Soc. B 289: 20220251 Compounded effects of climate change and habitat 279–295. (doi:10.1007/s00704-011-0481-1) comprehensive trait database of European and alteration shift patterns of butterfly diversity. Proc. 20. Pollard E. 1979 Population ecology and change in Maghreb butterflies, Papilionoidea. Sci. Data 7, Natl Acad. Sci. USA 107, 2088–2092. (doi:10.1073/ range of the white admiral butterfly Ladoga camilla 1–10. (doi:10.6084/m9.figshare.12998828) pnas.0909686107) L. in England. Ecol. Entomol. 4, 61–74. (doi:10. 34. Ninyerola M, Pons X, Roure JM. 2007 Objective air 7. Crossley MS. 2021 Recent climate change is creating 1111/j.1365-2311.1979.tb00561.x) temperature mapping for the Iberian Peninsula hotspots of butterfly increase and decline across 21. Van Nouhuys S, Lei G. 2004 Parasitoid–host using spatial interpolation and GIS. Int. J. Climatol. North America. Glob. Change Biol. 27, 2702–2714. metapopulation dynamics: the causes and 27, 1231–1242. (doi:10.1002/joc.1462) (doi:10.1111/gcb.15582) consequences of phenological asynchrony. J. Anim. 35. Moussus JP, Julliard R, Jiguet F. 2010 Featuring 8. Wilson RJ, Gutiérrez D, Gutiérrez J, Martínez D, Ecol. 73, 526–535. (doi:10.1111/j.0021-8790.2004. 10 phenological estimators using simulated data. Agudo R, Monserrat VJ. 2005 Changes to the 00827.x) Methods Ecol. Evol. 1, 140–150. (doi:10.1111/j. elevational limits and extent of species ranges 22. Altermatt F. 2010 Climatic warming increases 2041-210X.2010.00020.x) associated with climate change. Ecol. Lett. voltinism in European butterflies and moths. 36. Davis CC, Willis CG, Primack RB, Miller-Rushing AJ. 8, 1138–1146. (doi:10.1111/j.1461-0248.2005. Proc. R. Soc. B 277, 1281–1287. (doi:10.1098/rspb. 2010 The importance of phylogeny to the study of 00824.x) 2009.1910) phenological response to global climate change. 9. Pöyry J, Luoto M, Heikkinen RK, Kuussaari M, 23. Kerr NZ, Wepprich T, Grevstad FS, Dopman EB, Chew Phil. Tran. R. Soc. B 365, 3201–3213. (doi:10.1098/ Saarinen K. 2009 Species traits explain recent range FS, Crone EE. 2020 Developmental trap or rstb.2010.0130) shifts of Finnish butterflies. Glob. Change Biol. 15, demographic bonanza? Opposing consequences of 37. Wiemers M, Chazot N, Wheat CW, Schweiger O, 732–743. (doi:10.1111/j.1365-2486.2008.01789.x) earlier phenology in a changing climate for a Wahlberg N. 2020 A complete time-calibrated 10. Stefanescu C, Peñuelas J, Filella I. 2003 Effects of multivoltine butterfly. Glob. Change Biol. 26, multi-gene phylogeny of the European butterflies. climatic change on the phenology of butterflies in 2014–2027. (doi:10.1111/gcb.14959) ZooKeys 938, 97. (doi:zookeys.938.50878) the northwest Mediterranean Basin. Glob. Change 24. Kharouba HM, Vellend M. 2015 Flowering time of 38. R Core Team. 2021 R: a language and environment Biol. 9, 1494–1506. (doi:10.1046/j.1365-2486.2003. butterfly nectar food plants is more sensitive to for statistical computing. Vienna, Austria: R 00682.x) temperature than the timing of butterfly adult Foundation for Statistical Computing. https://www. 11. Diamond SE, Frame AM, Martin RA, Buckley LB. flight. J. Anim. Ecol. 84, 1311–1321. (doi: 10.1111/ R-project.org/. 2011 Species’ traits predict phenological responses 1365-2656.12373) 39. Wood SN. 2011 Fast stable restricted maximum to climate change in butterflies. Ecology 92, 25. Van Dyck H, Bonte D, Puls R, Gotthard K, Maes D. likelihood and marginal likelihood estimation of 1005–1012. (doi:10.1890/10-1594.1) 2015 The lost generation hypothesis: could climate semiparametric generalized linear models. J. R. Stat. 12. Macgregor CJ et al. 2019 Climate-induced change drive ectotherms into a developmental trap? Soc. B 73, 3–36. (doi:10.1111/j.1467-9868.2010. phenology shifts linked to range expansions in Oikos 124, 54–61. (doi:10.1111/oik.02066) 00749.x) species with multiple reproductive cycles per year. 26. Büntgen U, Piermattei A, Krusic P J, Esper J, 40. Bates D, Mächler M, Bolker B, Walker S. 2015 Fitting Nat. Commun. 10, 1–10. (doi:10.1038/s41467-019- Sparks T, Crivellaro A. 2022 Plants in the UK linear mixed-effects models using lme4. J. Stat. 12479-w) flower a month earlier under recent warming. Soft. 67, 1–48. (doi:10.18637/jss.v067.i01) 13. Hällfors MH, Pöyry J, Heliölä J, Kohonen I, Kuussaari Proc. R. Soc. B 289, 2021–2456. (doi: 10.1098/rspb. 41. Paradis E, Schliep K. 2019 ape 5.0: An environment M, Leinonen R, Schmucki R, Sihvonen P, 2021.2456) for modern phylogenetics and evolutionary analyses Saastamoinen M. 2021 Combining range and 27. Pollard E, Yates TJ. 1993 Monitoring butterflies for in R. Bioinformatics 35, 526–528. (doi:10.1093/ phenology shifts offers a winning strategy for boreal ecology and conservation: the British butterfly bioinformatics/bty633) Lepidoptera. Ecol. Lett. 24, 1619–1632. (doi:10. monitoring scheme. Berlin, Germany: Springer. 42. Barton K. 2020 MuMIn: multi-model inference. R 1111/ele.13774) 28. Schmucki R et al. 2016 A regionally informed package version 1.43.17. See https://CRAN.R-project. 14. Forister ML, Shapiro AM. 2003 Climatic trends and abundance index for supporting integrative analyses org/package=MuMIn. advancing spring flight of butterflies in lowland across butterfly monitoring schemes. J. Appl. Ecol. 43. Cohen JM, Lajeunesse MJ, Rohr JR. 2018 A global California. Glob. Change Biol. 9, 1130–1135. (doi:10. 53, 501–510. (doi:10.1111/1365-2664.12561) synthesis of animal phenological responses to 1046/j.1365-2486.2003.00643.x) 29. Melero Y, Stefanescu C, Pino J. 2016 General climate change. Nat. Clim. Change 8, 224–228. (doi: 15. Gutiérrez D, Wilson RJ. 2020 Intra-and interspecific declines in Mediterranean butterflies over the last 10.1038/s41558-018-0067-3) variation in the responses of insect phenology to two decades are modulated by species traits. Biol. 44. Lionello P, Scarascia L. 2018 The relation between climate. J. Anim. Ecol. 90, 248–259. (doi:10.1111/ Conserv. 201, 336–342. (doi:10.1016/j.biocon.2016. climate change in the Mediterranean region and 1365-2656.13348) 07.029) global warming. Reg. Environ. Change 18, 16. Parmesan C. 2007 Influences of species, latitudes 30. Zografou K, Swartz MT, Adamidis GC, Tilden VP, 1481–1493. (doi:10.1007/s10113-018-1290-1) and methodologies on estimates of phenological Mckinney EN, Sewall BJ. 2021 Species traits affect 45. Warren MS, Maes D, Van Swaay CA, Goffart P, Van response to global warming. Glob. Change Biol. phenological responses to climate change in a Dyck H, Bourn NA, Wynhoff I, Hoare D, Ellis S. 2021 The decline of butterflies in Europe: problems, 49. Michielini JP, Dopman EB, Crone EE. 2021 Changes change. Nat. Clim. Change 4, 127–131. (doi:10. 9 significance, and possible solutions. Proc. Natl Acad. in flight period predict trends in abundance of 1038/nclimate2086) royalsocietypublishing.org/journal/rspb Sci. USA 118. (doi:10.1073/pnas.2002551117) Massachusetts butterflies. Ecol. Lett. 24, 249–257. 53. Iler AM, Caradonna PJ, Forrest JR, Post E. 2021 46. Warren MS et al. 2001 Rapid responses of British (doi:10.1111/ele.13637) Demographic consequences of phenological shifts in butterflies to opposing forces of climate and habitat 50. Visser ME, Gienapp P. 2019 Evolutionary and response to climate change. Annu. Rev. Ecol. Evol. change. Nature 414, 65–69. (doi:10.1038/ demographic consequences of phenological Syst. 52, 221–245. (doi:10.1146/annurev-ecolsys- 35102054) mismatches. Nat. Ecol. Evol. 3, 879–885. (doi:10. 011921-032939) 47. Breed GA, Stichter S, Crone EE. 2013 Climate-driven 1038/s41559-019-0880-8) 54. Peñuelas J, Filella I, Comas P. 2002 Changed changes in northeastern US butterfly communities. 51. Spence AR, Tingley MW. 2020 The challenge plant and animal life cycles from 1952 to 2000 in the Nat. Clim. Change 3, 142–145. (doi:10.1038/ of novel abiotic conditions for species Mediterranean region. Glob. Change Biol. 8, 531–544. nclimate1663) undergoing climate-induced range shifts. (doi:10.1046/j.1365-2486.2002.00489.x) 48. Wepprich T, Adrion JR, Ries L, Wiedmann J, Haddad Ecography 43, 1571–1590. (doi:10.1111/ 55. Colom P, Ninyerola M, Pons X, Traveset A, NM. 2019 Butterfly abundance declines over 20 ecog.05170) Stefanescu C. 2022 Phenological sensitivity and years of systematic monitoring in Ohio, USA. PLoS 52. Mair L, Hill JK, Fox R, Botham M, Brereton T, seasonal variability explain climate-driven trends in Proc. R. Soc. B 289: 20220251 ONE 14, e0216270. (doi:10.1371/journal.pone. Thomas CD. 2014 Abundance changes and habitat Mediterranean butterflies. Figshare. (doi:10.6084/ 0216270) availability drive species’ responses to climate m9.figshare.c.5954002)