ARTICLE https://doi.org/10.1038/s42003-022-03088-3 OPEN Local adaptation to climate anomalies relates to species phylogeny Yolanda Melero 1,2 ✉, Luke C. Evans1, Mikko Kuussaari3, Reto Schmucki 4, Constantí Stefanescu5, David B. Roy 4 & Tom H. Oliver 1 Climatic anomalies are increasing in intensity and frequency due to rapid rates of global change, leading to increased extinction risk for many species. The impacts of anomalies are likely to vary between species due to different degrees of sensitivity and extents of local adaptation. Here, we used long-term butterfly monitoring data of 143 species across six 1234567890():,; European bioclimatic regions to show how species’ population dynamics have responded to local or globally-calculated climatic anomalies, and how species attributes mediate these responses. Contrary to expectations, degree of apparent local adaptation, estimated from the relative population sensitivity to local versus global anomalies, showed no associations with species mobility or reproductive rate but did contain a strong phylogenetic signal. The existence of phylogenetically-patterned local adaptation to climate has important implica- tions for forecasting species responses to current and future climatic conditions and for developing appropriate conservation practices. 1 School of Biological Sciences, University of Reading, Whiteknights, PO Box 217, Reading, Berkshire RG6 6AH, UK. 2 CREAF, E08193, Bellaterra, Cerdanyola del Vallès, Catalonia, Spain. 3 Finnish Environment Institute (SYKE), Biodiversity Centre, Latokartanonkaari 11, FI-00790 Helsinki, Finland. 4 UK Centre for Ecology & Hydrology, Wallingford, Oxfordshire OX10 8BB, UK. 5 Museu de Ciències Naturals de Granollers, Francesc Macià 51, Granollers, Catalonia 08402, Spain. ✉email: yolanda.melero@reading.ac.uk COMMUNICATIONS BIOLOGY | (2022)5:143 | https://doi.org/10.1038/s42003-022-03088-3 | www.nature.com/commsbio 1 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-022-03088-3 C limate change threatens many species with population species with less gene flow between populations and with higher decline and higher extinction risk, not only through the reproductive rates. Therefore, we expected increased adaptation effects of increases in mean temperature1–3 but also due to to local weather anomalies for (i) less mobile species, and (ii) the increased intensity and frequency of climatic anomalies (i.e., species with higher reproductive rates (assessed in terms of the variations from average climatic values4,5). For several species, number of generations per year). We found a wide range of climatic anomalies may have larger impacts than continuous responses across species, with some more sensitive to climatic increases in global mean temperature5,6. anomalies at the local site level, to others more sensitive to cli- Bioclimatic studies have reported population responses to matic anomalies calculated relative to a global mean across the micro-to-macro scaled climatic events7,8, but none have tested for whole dataset. We also found this interspecific variation was not contrasting population dynamics across the species distributional related to species mobility or reproductive rate, although it was range due to adaptations to the local biotic and abiotic conditions. phylogenetically patterned. Here, we hypothesize a continuum describing the degree of local adaptation for different species, from those showing relatively Results and discussion higher sensitivity to anomalies in the local-versus-global climatic Population change and climatic anomalies. To test our regime (“locally adapted” species) to those with higher sensitivity hypotheses, we used count data from 143 species collected across to global anomalies (“globally adapted”). This corresponds to the 172 sites within six European bioclimatic regions between 1999 known trade-offs between performance and tolerance, with and 2017 (167,288 site-year-species data points). For each site, populations of locally adapted species performing best at the butterfly counts were extracted from one of three European average local climatic conditions, but being sensitive to the local butterfly-monitoring schemes, all using a standardized transect- anomalies, hence showing population declines below and above based methodology (the “Pollard Walk”28). We used the weekly the optima; while globally adapted species will best perform counts collected across six bioclimatic regions to derive species’ around the average conditions occurring across their distribu- regional phenologies (flight period) and calculate series of annual tional range, and being sensitive more to the anomalies around abundance indices for each species and each site29. Specifically, this average (Fig. 1). Understanding the degree of local adaptation while accounting for density dependence (known to be important is essential for predicting species responses to climate change9,10. in butterfly populations30–32) we modeled how annual population To date, most research has focused on global or continental change per species varied in response to climatic anomalies at two species responses to climate change11,12, while local responses spatial scales: local (i.e., climatic anomalies in relation to the have received less attention13,14. Comprehension of species’ average across all years at that specific site) and global (i.e., cli- spatial adaptations at the appropriate scale is needed and likely to matic anomalies in relation to the average climate across all years lead to a better understanding of current species trends and and all sites where the species is found). improve predictions from bioclimatic models. Climatic anomalies were calculated for temperature, precipita- Previous research on the phenological responses of butterflies tion, and aridity occurring during species-specific phenological provides some evidence of local adaptation across species13, with periods (overwinter, pre-flight, flight, or post-flight period), both similar findings from single-species studies15. It is likely that for the year of population change t (except post-flight period) and species attributes contribute strongly to adaptive potential9,13 and the previous year t – 1, thereby accounting for possible direct and subsequent responses to global change16,17. Therefore, we expect delayed population responses. Through AIC-based model selec- that life-history traits mediate the degree of local adaptation to tion (Supplementary Data 1), we identified 86 species for which climatic anomalies. Experimentally designed research on insects population dynamics were significantly explained by climatic has shown that local adaptation to climate is mediated by both anomalies, of which half were affected by anomalies in physiological18,19 and morphological attributes (e.g., fur on body temperature (n = 48), either at the local or to the global scale, and melanin pigmentation20). However, it is unclear how these and 38 species responded to anomalies in precipitation or aridity. results translate to species population dynamics at broad spatial This result unifies findings from previous studies which have scales; partly due to the constraints of gathering long-term data related changes in populations of butterfly species individually to across the large distributional ranges of multiple species14. either temperature25, precipitation21, or aridity33,34. Overall, the In this study, we performed a comparative analysis using mul- importance of climatic anomalies at the most sensitive pheno- tiple species of butterflies as a study system because, like many other logical periods varied between the species, indicating idiosyncratic thermophilous organisms, they are highly sensitive to weather and interspecific variation in weather sensitivity2,5 (Supplementary strongly affected by climate change21–23, and show rapid responses Data 2). to changes24, minimizing (in comparison to plants and birds, for example) the demographic time lag to extinction25,26. Moreover, they show broad variation in life-history traits27, and have been Interspecific variation and the mediating role of species attri- subject to spatially- and temporally replicated standardized abun- butes. To test our first hypothesis, we quantified the extent to dance sampling across a large continental extent. which species were either locally or globally adapted by com- We first evaluate if species responses can be classified on our paring the variance in year-to-year population change explained hypothesized continuum of local-global adaptation. Based on a (R2) for each species using either local or global climatic simulation that demonstrates our hypothesis (Fig. 1), we expect anomalies, using the most appropriate climate variable for that that species adapted to local climatic conditions have population species (see above). The degree of apparent local adaptation for dynamics best explained by local rather than global climatic each species was calculated as the difference between the two R2 anomalies, i.e., they are more sensitive to anomalies occurring at values, ranging from −1 to 1, with higher values indicating the local spatial scale in terms of higher goodness of fit between greater relative sensitivity to local climatic conditions, and sug- interannual population change versus anomaly from average local gestive of local adaptation (Fig. 1). Despite some species showing climate conditions. Second, we test if local or global species’ a low degree of apparent local adaptation (<0.05), we observed a responses to climatic anomalies are mediated through key but- continuum of adaptation between species, from most locally terfly species attributes. adapted (e.g., Brenthis ino, Melitaea parthenoides) to most glob- We hypothesized that the relative sensitivity to local climate ally adapted (e.g., Cupido osiris, Apatura ilia; Fig. 2; Supple- anomalies (our proxy of local adaptation) would be higher for mentary Data 1 and Data 2). Species at the ends of the continuum 2 COMMUNICATIONS BIOLOGY | (2022)5:143 | https://doi.org/10.1038/s42003-022-03088-3 | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-022-03088-3 ARTICLE (a) (b) (c) 1 1 3 Population change Population change 0 0 Performance 2 -1 -1 -2 -2 1 -3 -3 0 -4 -4 0.00 0.25 0.50 0.75 1.00 -0.2 -0.1 0.0 0.1 0.2 -0.3 0.0 0.3 0.6 (d) (e) (f) 3 0 0 Population change Population change Performance -2 -2 2 -4 -4 1 -6 -6 0 -8 -8 0.00 0.25 0.50 0.75 1.00 -0.2 -0.1 0.0 0.1 0.2 -0.3 0.0 0.3 0.6 Scaled Climatic Variable Local Climatic Anomaly Global Climatic Anomaly Fig. 1 Simulations of the consequences of global and local adaptation on the population responses to local and global climatic anomalies. a, d show the performance of species at five sites with climatic means spanning across the range of the climatic variable. We expected locally adapted species to present multiple different performance curves representing distinct populations at sites distributed along the species’ distributional range, as shown in panel a. This expectation implies that population change will be more sensitive to local weather anomalies (simulation in b) than to weather anomalies calculated from all sites across the species’ distribution (simulation in c). In the case of global adaptation, performance is represented by a single curve through its entire range (d). Therefore, observed population change will be more sensitive to global weather anomalies calculated from all sites across the species’ distribution (f) than to the local site anomalies (e). Performance curves were based on a Briere type I function74, which is a simple function that matches empirical data on thermal performance75. We included a fixed area under the curve as consistent with expectations of specialist/generalist trade-offs76. Beyond this hypothetical example, in practice, the mean and variance of curves may vary across species; for example, some species may have broader climatic tolerances than others, i.e., the curves in panels b or f would be broader and shallower, meaning there is a greater range of temperatures in which population growth can remain positive. Broader tolerances may be driven in part by phenotypic plasticity, i.e., gene-by-environment interactions (for example, oviposition microsite preferences varying between locations depending on the local macroclimate). This phenotypic plasticity may be exhibited across the entire range, or it might only occur in certain areas, i.e., there might conceivably be a local adaptation of the phenotypic plasticity. Alternatively, local adaptation can occur in fixed traits, such as lighter-colored insects in warmer areas77. Both of these evolutionary adaptations produce patterns akin to that in panel a, whereby optimum of a thermal performance differs amongst the populations of species so that they perform best in their “home” conditions. To generate weather across the range, we standardized an observed 19-year time series of global yearly temperatures (min = 0, max = 1, mean = 0.5) and then shifted the values of each year to predict mean expectation at local sites across the range, a local value for each site and year was then sampled with Gaussian noise. The performance was subsequently used as input into a discrete logistic growth model (Nt+1 = RNt(1 – Nt/K)) as proportional to R the intrinsic growth rate. Each population was seeded with a small number of individuals and was allowed to recover by immigration should the population size go to zero. A time series of population change for each of the sites was collected from the simulation (ΔN after initialization and immigration was excluded). Models for population change were then fitted using local and global anomalies and are shown in (b, c, e, f). Colors in (a, d) indicate location in the distributional (e.g., blue to red, cold to hot extremes, respectively). Colors in the rest of the panels indicate spatial scale (blue—local climate anomalies; orange—global climate anomalies), circles indicate populations (i.e., from distinct sites), lines show predicted trend with 95% confident intervals. responded much more strongly to anomalies at a specific scale their phylogenetic relationships (see “Methods”). We found a (local or global), producing unimodal “n-shaped” responses, strong phylogenetic signal, explaining 84% of the variation which are expected given sampling across their climatic niche (Pagel’s lambda = 0.84; P value < 0.0001; Fig. 4), suggesting (illustrative examples provided for Brenthis ino and Cupido osiris phylogenetic patterning in the degree of local adaptation across Fig. 3; Supplementary Figs. 1–86). species (Fig. 1). However, we found no effect of either mobility For the second hypothesis, we tested if relative sensitivity to or reproductive rate (Supplementary Fig. 87 and Supplementary local-versus-global climate anomalies was mediated by species’ Table 1); although mobility and reproductive rate had a mobility35 and reproductive rate (number of generations per significant phylogenetic signal, this did not explain the degree year; i.e., voltinism, with species categorized as univoltine versus of species local adaptation beyond that expected from multivoltine36). We also accounted for the potential effect of taxonomic relatedness (Supplementary Fig. 88 and COMMUNICATIONS BIOLOGY | (2022)5:143 | https://doi.org/10.1038/s42003-022-03088-3 | www.nature.com/commsbio 3 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-022-03088-3 Fig. 2 Degree of apparent local adaptation across European butterfly species, determined by relative sensitivity of population dynamics to local- versus-global climatic anomalies (n = 86). Degree of local adaptation is calculated as the difference between the R2 of the model including the climatic anomaly at the local scale (site-specific) and the model including that correspondent climatic anomaly at the global scale (all sites), with positive values (blue points) indicating greater local adaption and orange points indicating greater global adaptation. Letters indicate the climatic variable affecting the species dynamics: T temperature, P precipitation, A aridity. Supplementary Table 2). Given these results, it was possible that reproductive rates. Note, that there was no effect of GST when species overall sensitivity to climatic anomalies (either local or testing it as a proxy of genetic differentiation as an alternative to global) or lack of sensitivity (e.g., species for which density mobility (Supplementary Table 3). Phylogenetic signals have also dependence was a main factor), also had a phylogenetic signal. been seen in distributional shifts of plant communities40 and Hence, we performed a post hoc analysis on the 143 species and some butterfly species41; while mixed results have been found in found that phylogeny explained 88% of the species presence or relation to mobility and reproductive rate in range shifts or lack of sensitivity to climatic anomalies (Supplementary Fig. 89 distributional trends42,43. Specific phylogenetic patterns of local and Supplementary Table 2). The strong phylogenetic signal of adaptation might also result from the pressures of the external both the degree of local adaptation (relative sensitivity to local- conditions affecting the capacity to evolve and adapt. For versus-global climate anomalies) and the overall total sensitivity example, local adaptation might be more likely in certain regions to climatic anomalies, provides robustness to discard potential or habitat types, with phylogenetic patterning in the associations stochasticity in the results. of species to these conditions44,45. Globally adapted species could result from low local selective pressures, although we found no association between the overall The role of species phylogeny mediating local adaptation. Our variance in population dynamics explained by climate with results suggest a dominant effect of evolutionary constraints (i.e., species degree of local-versus-global adaptation (P value = 0.54; unmeasured traits with a strong phylogenetic signal) on the Supplementary Fig. 90). Alternatively, global adaptation might capacity for local adaptation. However, they provide no support result from a lack of genetic variation across the populations46, or for our expectation of low gene flow and high reproductive rates could be linked to high gene flow that is not reflected, or leading to increased adaptive potential. remained undetected, in our measure of species dispersal An apparent degree of local adaptation shared by closely ability47. Extensive gene flow between populations and weak related species might be expected. Greater sensitivity to local selection selective pressures may constrain local adaptation48, climate anomalies implies spread and survival at sites where though recent studies indicate that this is not always the case for adaptation to local conditions is selected; this could be facilitated gene flow49,50. Finally, global adaptation may result from by the filtering of individuals from pre-existing genetic diversity pleiotropy, i.e., single genes affecting more than one phenotypic or through phenotypic plasticity (e.g., behavioral, morphological trait, limiting the degree of adaptive change due to selective trade- or physiological attributes) favouring survival and eventual offs across the attributes for some species51,52. evolutionary adaptation. Closely related species share similar genetic and phenotypic attributes, although we did not find mediation of local adaptation through either mobility or Implications for ecological forecasting and conservation. It is voltinism. This is possibly because European butterfly species important to contextualize the intra- and interspecific variation of have resided in their current location for long periods (many population fluctuations of butterfly species when facing climatic butterflies expanded across Europe during their post-glacial anomalies because they have implications for ecological fore- colonization37,38), even creating refugia in Southern areas39, casting and conservation under climate change9,53. Local adap- resulting in a long-term established genetic differentiation tations are likely to enhance stability for some populations under currently poorly affected by mobility or reproductive rates, and perturbations. For example, assuming similar niche breadth, hence with ample time to adapt to local conditions regardless of locally adapted species at their retracting warm distributional 4 COMMUNICATIONS BIOLOGY | (2022)5:143 | https://doi.org/10.1038/s42003-022-03088-3 | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-022-03088-3 ARTICLE 3 3 (a) R2 = 0.15* (b) R2 = 0.08 2 2 1 1 Population Change 0 0 −1 −1 −2 −2 −3 −3 −6 −4 −2 0 2 4 −5 −4 −3 −2 −1 0 3 3 (c) R2 = 0.24 (d) R2 = 0.47* 2 2 1 1 Population Change 0 0 −1 −1 −2 −2 −3 −3 −3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3 Local Climatic Anomaly Global Climatic Anomaly Fig. 3 Population change in relation to local and global climatic anomalies. a, b show local and global responses, respectively, for Brenthis ino, a species best adapted to local climatic anomalies in temperature during the overwintering period of the previous year of their adult stage (t − 1). c, d show local and global adaptation for Cupido osiris, a species whose population dynamics are best explained by global climatic anomalies in aridity during the pre-flight period of the year of their adult stage (t), indicating a lack of local adaptation. Colors indicate spatial scale (blue, local; orange, global), circles indicate raw data, lines show predicted trend with 95% confident intervals, asterisks indicate the best model for each species. edges may remain within their fundamental niche space for conditions over just one or two generations, potentially mediated longer than expected, compared with if they were globally by epigenetic changes58. Nonetheless, the relatively strong adapted. This difference in species capacity to cope with envir- phylogenetic signal found here is contingent upon a few species onmental change could be of particular importance to the con- with large leverage. Removing species based on data availability, servation of globally adapted species if populations at warm or based on visualization of outliers in the tree, leads to a lower marginal edges contain the majority of the species genetic phylogenetic signal (Supplementary Fig. 91). We also do not diversity54. In contrast, locally adapted populations at the colder discard the role of other factors such as habitat availability and leading edges of their ranges, expected to be strongholds of per- connectivity in influencing butterfly species dynamics as well as sistence under climate warming, may be more vulnerable to local potential interactions with the effect of the climatic anomalies, heat extremes than if they were globally adapted. We suggest but their role was out of the scope of this study. future research could seek to quantify if or how local adaptation Despite raising interesting questions to pursue further, our interacts with the specific location of populations within their overall results hold important implications for forecasting species geographic range to influence population changes, also taking responses to climate change. Namely, that local adaptation is account of how climatic sensitivity within species can vary across important to factor into predictive climate models for a substantial their geographical range55,56. subset of species, and this is likely to also influence advice We acknowledge some limitations of the study, with some regarding appropriate conservation practices in different locations. species (such as Vanessa cardui, Supplementary Fig. 86) showing unexpected strong signals of local adaptation despite multi- generational trans-continental migration of this species57. This Methods might be due to spurious results but could be genuine and due to Data collection. Butterfly count data were extracted from three long-term Butterfly factors that we cannot elucidate within this study but are worth Monitoring Schemes carried out in Finland, Spain, and the UK, covering six out of the ten bioclimatic regions across Europe and three main biomes (Supplementary researching in the future; for example, the potential existence of Figs. 92 and 93). The schemes consist of a network of sites where volunteers several races with different environmental tolerances and perform visual counts of butterflies along transects (i.e., sites) following the stan- dispersal patterns across Europe, or rapid adaptation to local dardized “Pollard Walk” methodology28. Counts are conducted weekly during the COMMUNICATIONS BIOLOGY | (2022)5:143 | https://doi.org/10.1038/s42003-022-03088-3 | www.nature.com/commsbio 5 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-022-03088-3 Fig. 4 Phylogenetic comparative tree for the degree of local adaptation for the 86 species significantly affected by climatic anomalies. Degree of local adaptation was set as a continuous distribution from −1 to 1. Negative (red) values of the degree of local adaptation relate to species adapted to climatic anomalies at the global scale, positive values (blue) the local scale. butterfly flight season, which varies depending on the climatic zone within the used this annual species index of abundance (Ni) to calculate the interannual range of March to the end of September. The species nomenclature is defined change in abundance as a proxy of the annual population growth rate. The rate following Wiemers et al.59. was calculated per species population (i.e., per site) between two consecutive For each site, year and species, we calculated an annual abundance index of years as log(Nt+1/Nt). adult butterflies that accounted for missing counts and spatial and temporal The duration of the schemes (i.e., starting year) and the number of surveyed variation in species phenology29,60. We derived the annual species phenology sites varied between countries: Finland (1999, nsites = 107), Spain (1994, from the observed weekly counts, fitting generalized additive models (GAM) for nsites = 130) and the UK. (1976, nsites = 2128). Thus, we used Finland (the limiting each bioclimatic region29. We then used the regional flight curve derived from country) to set the range of study years (1999–2017) and the number of sites per the GAM as an offset in a generalized linear model to predict and impute values country. To have sufficient interannual change data per site, we retained only sites for missing counts. The complete series of weekly counts, including the imputed having at least ten years of interannual abundance change data between 1999 and values, was then used to calculate an annual abundance index per species, site, 2017. This rule reduced our sample to 53 sites for Finland, 59 for Spain and 701 for and year. This method corrects for missing week counts, using the information the UK. From the 701 UK sites, we used a random subset of 60 sites to balance the derived from the regional phenology and produce indices that are less biased sample size between countries (Supplementary Figs. 94–99). These sites were than indices derived from linear interpolations29. However, to assure consistent selected randomly, apart from a condition maintaining at least 20 km separation estimates, we excluded indices of abundance >50% of missing observations. We between the selected sites. 6 COMMUNICATIONS BIOLOGY | (2022)5:143 | https://doi.org/10.1038/s42003-022-03088-3 | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-022-03088-3 ARTICLE As the effect of climate on population growth rates varies depending on the life variable, spatial scale (local or global), and time (t or t − 1) and phenological stages of the organisms61,62, we defined four life stages periods relevant for period, that most affected the species populations dynamics. Model selection was butterfly species55: pre-flight period, flight period, post-flight period, and based on AIC value (Supplementary Data 2). Species with the null model as the overwintering. Pre-, post-, and flight periods vary across species, years, and best model were discarded from all further analyses (i.e., AIC null model—AIC bioclimatic zone29. Thus, we set these periods per species, year and bioclimatic alternative model <2; n = 56). We did this on the basis that in these cases our zone62, except the post-flight period of year t, which we discarded because it could climatic variables explained almost nothing of the species population dynamics, i.e., have no possible effect on the population dynamics of the species in that year the variance explained by models with climatic variables was close to null. We also (adults in the count of that year had already perished). These periods were defined discarded one species (Polyommatus ripartii) as the ΔAIC between models using the annual flight curve distribution from the relative abundances per species containing different climatic anomalies, different periods, times and scale differed and zone (i.e., using the GAM fitted annual species phenology). The flight period <2; such as it was not possible to select the best model. We retained any species was set as the dates between the 10th and 90th percentiles of the species with the best models, including both the global and local anomalies of a specific distribution for each year and each bioclimatic zone. The pre-flight period was then climatic variable, at a specific period and time (n = 11). We did so because we were set from February to the 10th percentile of the species flight period, and the post- interested in the degree of local adaptation (see below); therefore both scales were flight period from the 90th percentile to the end of October. For multivoltine eventually used. Model selection led to 86 species for which the best models species (i.e., species with more than one reproduction per year), flight period was included climatic anomalies at the local and the global scale, all with AIC lower set between the 10th and 90th percentiles of the entire species flight-period than that of the null model. Following a conservative approach, we initially distribution, independently of the generation. The overwintering period was set as discarded nine out of these 86 species because they had less than ten years of fixed from November to January (of year t – 1 to year t), equally for all species and interannual abundance change data (Supplementary Data 1 and Data 2). all sites. Although overwintering could be defined over different time periods, we Discarding them did not influence the results related to attributes (see next fixed it to minimize potential overlaps between it and its adjacent time periods for section), but implied that some species were alone in their clade (e.g., Apatura ilia) some species in some regions55. or completely separated from the other groups (Papilio machaon). Therefore, we We calculated climatic anomalies for temperature, precipitation, and aridity (as decided to add them back to have a more complete and less biased phylogenetic a combination of temperature and precipitation) because previous studies have tree that encompassed all 86 species. demonstrated that butterflies, as also many other organisms, can be affected by each of these variables21,33,34,63. Climatic data were obtained from the European The role of phylogeny and species attributes mediating the effect of climatic Climatic and Assessment Dataset project (ECAD)64,65. The data consisted of daily anomalies in species population dynamics. For each species, we extracted the temperature (in degrees Celsius) and precipitation (in millimeters) from a gridded conditional R2, i.e., the variance explained by the fixed effects without accounting climatic dataset constructed from extrapolations of observations from a series of for the variance due to random effects of the best fitting model and that of the meteorological stations throughout Europe. We extracted the daily temperature model with the same climate variable but measured at the alternative spatial scale and precipitation at 1° scale for each site for the period 1999–2017; and calculated (local or global), keeping all other specifications (i.e., life-history period and year) local and global climatic anomalies for temperature and precipitation per species fixed. We then used these variances to calculate a proxy of the degree of local and phenological period (overwintering, pre-, post-, and flight period) as follows: adaptation per species by subtracting R2local model – R2global model. This proxy was Wsi; local ¼ Wsi  Ws ð1Þ subsequently used as the response variable to evaluate how local adaptation relates to mobility, voltinism, and phylogeny. To do so, we used a phylogenetic generalized Wsi; global ¼ Wsi  W ð2Þ linear model67, with the species degree of local adaptation fitted to a normal error distribution, and mobility and voltinism as explanatory variables while controlling where Wsi is the mean temperature or the total rainfall, depending on the climatic for the phylogenetic signal of the degree of local adaptation calculated using Pagel’s variable being analyzed per site (s) and year (i), Ws is the mean temperature or the lambda statistic68. total rainfall per site for all studied years, and W is the mean temperature or total Mobility and voltinism attribute values were retrieved from the literature35,36,69, rainfall for all studied sites and all years. Potential interactions between temperature with voltinism categorized as either univoltine or multivoltine species (i.e., strictly and rainfall were calculated as the standardized aridity index (Eqs. (3) and (4))33,66: one generation per year, and two or more generations per year respectively).  Species that have a rare, or occasional, second generation were set as univoltine. SAI si;local ¼ ððPsi  Ps Þ=sdPs Þ  0:5ðT si  T s Þ=sdT s ð3Þ The butterfly phylogenetic tree was extracted from the most recently published phylogenetic tree of European butterflies37. Alternative categorizations of voltinism SAI si;global ¼ ððPsi  PÞ=sdPÞÞ  0:5ðT si  TÞ=sdTÞ ð4Þ were also tested with species categorized as univoltine, putative multivoltine or where SAI stands for Standardized Aridity Index, T for mean temperature, P for strict multivoltine, and categorized as per their know number maximum total precipitation, sd for standard deviation, with subindexes as above. generations (1, 1.5, 2, or more than 2); we also included potential variability of voltinism set as the maximum and the minimum number of generation shown by a given species across Europe. Duration of flight period (in average months) and Modeling population changes in relation to climatic anomalies. We fitted three genetic differentiation (GST)37 were also tested as alternative proxies of genetic models to test whether species population growth rates were locally or globally differentiation and genetic flow (mobility). None of these alternatives differed adapted: a null model accounting only for density dependence to remove species qualitatively from the results of our initial model (Supplementary Table 3). not affected by the climatic anomalies when testing our hypotheses (Eq. (5)); a local A post hoc analysis for the phylogenetic signal of the species sensitivity to model accounting for density dependence, the local anomaly term of the climatic climatic anomalies was done using D statistic70; i.e., whether climate was important variable being analyzed, and its quadratic term (Eq. (6)); and a global model with in explaining the species population dynamics (best model included climatic density dependence, the global anomaly of the climatic variable, and its quadratic anomalies, N = 86) or whether density dependence best explained the population term (Eq. (7)). Models were fitted independently per each species, with population dynamics (the null model was the best model). growth rate as the response variable expressed as log(Nit/Nit − 1). Site was set as a We conducted all analyses in R 3.6.1 (R Core Team, 2019), using the package random intercept for all species, except for those with a low number of sites lme471 to fit our models and the package MuMin72 to estimate the marginal and (n = 23 species), due to lack of model convergence when adding random effect; in conditional variances. Phylogenetic analyses were done using the phylosig and pgls these cases, only fixed effects were added (Supplementary Data 2). Errors were functions in R package phytools73. assumed to be normally distributed, and model diagnostics were used to test the conformation of the data to the model assumptions. logðN it =N it1 Þ ¼ N it1 þ ε ð5Þ Data availability The count and abundance butterfly data that support the findings of this study are available from the European Butterfly Monitor Scheme via a signed license agreement logðN it =N it1 Þ ¼ N it1 þ Wsi;localt þ W2 si;localt þ ε ð6Þ (https://butterfly-monitoring.net/). Climatic data are available via ECAD website (https:// www.ecad.eu/). logðN it =N it1 Þ ¼ N it1 þ Wsi;globalt þ W2 si;globalt þ ε ð7Þ Where Ni is the annual abundance index at the site ith and time t or t – 1, Wit is the climatic anomaly variable site i at time t either at the local or the global scale. Received: 21 January 2021; Accepted: 28 January 2022; Models including climate were fitted per species with all combinations of climatic variable W (temperature, precipitation, and aridity), local or the global scale, time t or t − 1 (Wit and Wit-1, the latter not shown in the equations for simplicity), and all phenological periods (winter, pre-, post-, or flight period; except for the post-flight period at time t; Eqs. (6) and (7)): This procedure resulted in a set of 43 models per species: the reduced model References and 42 models consisting of each of the three climatic variables, per two spatial 1. Verdura, J. et al. Biodiversity loss in a Mediterranean ecosystem due to an scales (local and global), per two time periods (t and t − 1), and per four time extreme warming event unveils the role of an engineering gorgonian species. periods (excepting the post-flight period at time t). From within these 43 models, Sci. Rep. 9, 1–11 (2019). we retained the best fitting model for each species for identifying the climatic COMMUNICATIONS BIOLOGY | (2022)5:143 | https://doi.org/10.1038/s42003-022-03088-3 | www.nature.com/commsbio 7 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-022-03088-3 2. Pandori, L. L. M. & Sorte, C. J. B. The weakest link: sensitivity to climate 33. Oliver, T. H. et al. Interacting effects of climate change and habitat extremes across life stages of marine invertebrates. Oikos 128, 621–629 (2019). fragmentation on drought-sensitive butterflies. Nat. Clim. Chang. 5, 941–946 3. Palmer, G. et al. Climate change, climatic variation and extreme biological (2015). responses. Philos. Trans. R. Soc. B Biol. Sci. 372, 20160144 (2017). 34. Stefanescu, C., Carnicer, J. & Peñuelas, J. Determinants of species richness in 4. Altwegg, R., Visser, V., Bailey, L. D. & Erni, B. Learning from single extreme generalist and specialist Mediterranean butterflies: the negative synergistic events. Philos. Trans. R. Soc. B Biol. Sci. 372, 20160141 (2017). forces of climate and habitat change. Ecography 34, 353–363 (2011). 5. McDermott Long, O. et al. Sensitivity of UK butterflies to local climatic 35. Essens, T., van Langevelde, F., Vos, R. A., Van Swaay, C. A. M. & extremes: which life stages are most at risk? J. Anim. Ecol. 86, 108–116 (2017). WallisDeVries, M. F. Ecological determinants of butterfly vulnerability across 6. Jentsch, A., Kreyling, J. & Beierkuhnlein, C. A new generation of climate‐ the European continent. J. Insect Conserv. 21, 439–450 (2017). change experiments: events, not trends. Front. Ecol. Environ. 5, 365–374 36. Tolman, T. & Lewington, R. Butterflies of Europe (Harper Collins, 2008). (2007). 37. Dapporto, L. et al. Integrating three comprehensive data sets shows that 7. Suggitt, A. J. et al. Habitat associations of species show consistent but weak mitochondrial DNA variation is linked to species traits and paleogeographic responses to climate. Biol. Lett. 8, 590–593 (2012). events in European butterflies. Mol. Ecol. Resour. 19, 1623–1636 (2019). 8. Trisos, C. H., Merow, C. & Pigot, A. L. The projected timing of abrupt 38. Hewitt, G. M. Post-glacial re-colonization of European biota. Biol. J. Linn. Soc. ecological disruption from climate change. Nature 580, 496–501 (2020). Lond. 68, 87–112 (2008). 9. Valladares, F. et al. The effects of phenotypic plasticity and local adaptation on 39. Dincă, V. et al. High resolution DNA barcode library for European butterflies forecasts of species range shifts under climate change. Ecol. Lett. 17, reveals continental patterns of mitochondrial genetic diversity. Commun. Biol. 1351–1364 (2014). 4, 1–11 (2021). 10. Bush, A. et al. Incorporating evolutionary adaptation in species distribution 40. Fei, S. et al. Divergence of species responses to climate change. Sci. Adv. 3, modelling reduces projected vulnerability to climate change. Ecol. Lett. 19, e1603055 (2017). 1468–1478 (2016). 41. Macgregor, C. J. et al. Climate-induced phenology shifts linked to range 11. Stephens, P. A. et al. Consistent response of bird populations to climate expansions in species with multiple reproductive cycles per year. Nat. change on two continents. Science 352, 84–87 (2016). Commun. 10, 1–10 (2019). 12. Kerr, J. T. et al. Climate change impacts on bumblebees converge across 42. Dapporto, L. & Dennis, R. L. H. The generalist–specialist continuum: testing continents. Science 349, 177–180 (2015). predictions for distribution and trends in British butterflies. Biol. Conserv. 13. Roy, D. B. et al. Similarities in butterfly emergence dates among 157, 229–236 (2013). populations suggest local adaptation to climate. Glob. Chang. Biol. 21, 43. MacLean, S. A. & Beissinger, S. R. Species’ traits as predictors of range shifts 3313–3322 (2015). under contemporary climate change: a review and meta-analysis. Glob. Chang. 14. Titeux, N. et al. The need for large-scale distribution data to estimate regional Biol. 23, 4094–4105 (2017). changes in species richness under future climate change. Divers. Distrib. 23, 44. Morlon, H. et al. Spatial patterns of phylogenetic diversity. Ecol. Lett. 14, 1393–1407 (2017). 141–149 (2011). 15. Haeler, E., Fiedler, K. & Grill, A. What prolongs a butterfly’s life?: trade-offs 45. Kraft, N. J. B. et al. Community assembly, coexistence and the environmental between dormancy, fecundity and body size. PLoS One 9, e111955 (2014). filtering metaphor. Funct. Ecol. 29, 592–599 (2015). 16. Gonzalez-Suarez, M., Gomez, A. & Revilla, E. Which intrinsic traits predict 46. Razgour, O. et al. Considering adaptive genetic variation in climate change vulnerability to extinction depends on the actual threatening processes. vulnerability assessment reduces species range loss projections. Proc. Natl Ecosphere 4, 1–16 (2013). Acad. Sci. USA 116, 10418–10423 (2019). 17. Pacifici, M. et al. Species’ traits influenced their response to recent climate 47. Vanden Broeck, A. et al. Gene flow and effective population sizes of the change. Nat. Clim. Chang. 7, 205–208 (2017). butterfly Maculinea alcon in a highly fragmented, anthropogenic landscape. 18. Kingsolver, J. G. & Watt, W. B. Thermoregulatory strategies in Colias Biol. Conserv. 209, 89–97 (2017). butterflies: thermal stress and the limits to adaptation in temporally varying 48. Haldane, J. B. S. Theoretical genetics of autopolyploids. J. Genet. 22, 359–372 environments (Colorado). Am. Nat. 121, 32–55 (1983). (1930). 19. MacLean, H. J., Higgins, J. K., Buckley, L. B. & Kingsolver, J. G. Morphological 49. Tigano, A. & Friesen, V. L. Genomics of local adaptation with gene flow. Mol. and physiological determinants of local adaptation to climate in Rocky Ecol. 25, 2144–2164 (2016). Mountain butterflies. Conserv. Physiol. 4, 1 (2016). 50. Pfeifer, S. P. et al. The evolutionary history of Nebraska deer mice: local 20. Kingsolver, J. G. & Wiernasz, D. C. Seasonal polyphenism in wing-melanin adaptation in the face of strong gene flow. Mol. Biol. Evol. 35, 792–806 (2018). pattern and thermoregulatory adaptation in Pieris butterflies. Am. Nat. 137, 51. Reusch, T. B. H. & Wood, T. E. Molecular ecology of global change. Mol. Ecol. 816–830 (1991). 16, 3973–3992 (2007). 21. Herrando, S. et al. Contrasting impacts of precipitation on Mediterranean 52. DeLong, J. P. & Gibert, J. P. Gillespie eco-evolutionary models (GEMs) reveal birds and butterflies. Sci. Rep. 9, 1–7 (2019). the role of heritable trait variation in eco-evolutionary dynamics. Ecol. Evol. 6, 22. Thomas, J. A. Monitoring change in the abundance and distribution of insects 935–945 (2016). using butterflies and other indicator groups. Philos. Trans. R. Soc. Lond. B. 53. Atkins, K. E. & Travis, J. M. J. Local adaptation and the evolution of species’ Biol. Sci. 360, 339–357 (2005). ranges under climate change. J. Theor. Biol. 266, 449–457 (2010). 23. Roy, D. B., Rothery, P., Moss, D., Pollard, E. & Thomas, J. A. Butterfly 54. Hampe, A. & Petit, R. J. Conserving biodiversity under climate change: the numbers and weather: predicting historical trends in abundance and the rear edge matters. Ecol. Lett. 8, 461–467 (2005). future effects of climate change. J. Anim. Ecol. 70, 201–217 (2008). 55. Mills, S. C. et al. European butterfly populations vary in sensitivity to weather 24. Pöyry, J., Luoto, M., Heikkinen, R. K., Kuussaari, M. & Saarinen, K. Species across their geographical ranges. Glob. Ecol. Biogeogr. 26, 1374–1385 (2017). traits explain recent range shifts of Finnish butterflies. Glob. Chang. Biol. 15, 56. Van Dyck, H., Bonte, D., Puls, R., Gotthard, K. & Maes, D. The lost generation 732–743 (2009). hypothesis: could climate change drive ectotherms into a developmental trap? 25. Devictor, V. et al. Differences in the climatic debts of birds and butterflies at a Oikos 124, 54–61 (2015). continental scale. Nat. Clim. Chang. 2, 121–124 (2012). 57. Hu, G. et al. Environmental drivers of annual population fluctuations in a 26. Krauss, J. et al. Habitat fragmentation causes immediate and time-delayed trans-Saharan insect migrant. Proc. Natl Acad. Sci. USA 118, 2102762118 biodiversity loss at different trophic levels. Ecol. Lett. 13, 597–605 (2010). (2021). 27. Eskildsen, A. et al. Ecological specialization matters: long-term trends in 58. Merlin, C. & Liedvogel, M. The genetics and epigenetics of animal migration butterfly species richness and assemblage composition depend on multiple and orientation: birds, butterflies and beyond. J. Exp. Biol. 222, jeb191890 functional traits. Divers. Distrib. 21, 792–802 (2015). (2019). 28. Pollard, E. A method for assessing changes in the abundance of butterflies. 59. Wiemers, M. et al. An updated checklist of the European butterflies Biol. Conserv. 12, 115–134 (1977). (Lepidoptera, Papilionoideae). Zookeys 2018, 9–45 (2018). 29. Schmucki, R. et al. A regionally informed abundance index for supporting 60. Dennis, E. B., Freeman, S. N., Brereton, T. & Roy, D. B. Indexing butterfly integrative analyses across butterfly monitoring schemes. J. Appl. Ecol. 53, abundance whilst accounting for missing counts and variability in seasonal 501–510 (2016). pattern. Methods Ecol. Evol. 4, 637–645 (2013). 30. Pollard, E., Lakhani, K. H. & Rothery, P. The detection of density-dependence 61. Radchuk, V., Turlure, C. & Schtickzelle, N. Each life stage matters: the from a series of annual censuses. Ecology 68, 2046–2055 (1987). importance of assessing the response to climate change over the complete life 31. Dooley, C. A., Bonsall, M. B., Brereton, T. & Oliver, T. Spatial variation cycle in butterflies. J. Anim. Ecol. 82, 275–285 (2013). in the magnitude and functional form of density-dependent processes on 62. Metzger, M. J. et al. A high-resolution bioclimate map of the world: a unifying the large skipper butterfly Ochlodes sylvanus. Ecol. Entomol. 38, 608–616 framework for global biodiversity research and monitoring. Glob. Ecol. (2013). Biogeogr. 22, 630–638 (2013). 32. Rothery, P., Newton, I., Dale, L. & Wesolowski, T. Testing for density 63. Carnicer, J. et al. A unified framework for diversity gradients: the adaptive trait dependence allowing for weather effects. Oecologia 112, 518–523 (1997). continuum. Glob. Ecol. Biogeogr. 22, 6–18 (2013). 8 COMMUNICATIONS BIOLOGY | (2022)5:143 | https://doi.org/10.1038/s42003-022-03088-3 | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-022-03088-3 ARTICLE 64. Klok, E. J. & Klein Tank, A. M. G. Updated and extended European dataset of Author contributions daily climate observations. Int. J. Climatol. 29, 1182–1191 (2009). L.E., M.K., D.R., R.S., C.S., T.O., and Y.M. contextualized the study and its 65. Haylock, M. R. et al. A European daily high-resolution gridded data set of hypotheses. Y.M. led the paper and performed the analyses, with L.E. and surface temperature and precipitation for 1950-2006. J. Geophys. Res. Atmos. supervised by T.O. L.E., M.K., D.R., R.S., C.S., T.O., and Y.M. all contributed to 113, D20119 (2008). writing the manuscript. 66. Marsh, T. J. The UK drought of 2003: a hydrological review. Weather 59, 224–230 (2004). 67. Voyer, A. G. & Garamszegi, L. Z. An introduction to phylogenetic path Competing interests The authors declare no competing interests. analysis. in Modern Phylogenetic Comparative Methods and their Application in Evolutionary Biology (eds Garamszegi, L. Z. & Mundry, R.) 201–229 (Springer Berlin Heidelberg, 2014). Additional information 68. Pagel, M. Inferring the historical patterns of biological evolution. Nature 401, Supplementary information The online version contains supplementary material 877–884 (1999). available at https://doi.org/10.1038/s42003-022-03088-3. 69. Pöyry, J. et al. The effects of soil eutrophication propagate to higher trophic levels. Glob. Ecol. Biogeogr. 26, 18–30 (2017). Correspondence and requests for materials should be addressed to Yolanda Melero. 70. Münkemüller, T. et al. How to measure and test phylogenetic signal. Methods Ecol. Evol. 3, 743–756 (2012). Peer review information Communications Biology thanks Leonardo Dapporto and the 71. Bates, D., Maechler, M., Bolker, B. & Walker, S. Fitting Linear Mixed-Effects other, anonymous, reviewers for their contribution to the peer review of this work. Models Using lme4. J. Stat. Softw. 67, 1–48 (2015). Primary Handling Editors: Michelle Lawing and Caitlin Karniski. Peer reviewer reports 72. Bartoń, K. MuMIn: Multi-model inference. R package version 1.10.5. (2014). are available. 73. Revell, L. J. phytools: An R package for phylogenetic comparative biology (and other things). MEE. 3, 217–223 (2012). Reprints and permission information is available at http://www.nature.com/reprints 74. Briere, J. F., Pracros, P., Le Roux, A. Y. & Pierre, J. S. A novel rate model of temperature-dependent development for arthropods. Environ. Entomol. 28, Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in 22–29 (1999). published maps and institutional affiliations. 75. Shi, P. & Ge, F. A comparison of different thermal performance functions describing temperature-dependent development rates. J. Therm. Biol. 35, 225–231 (2010). Open Access This article is licensed under a Creative Commons 76. Angilletta, M. J., Wilson, R. S., Navas, C. A. & James, R. S. Tradeoffs and the Attribution 4.0 International License, which permits use, sharing, evolution of thermal reaction norms. Trends Ecol. Evol. 18, 234–240 (2003). adaptation, distribution and reproduction in any medium or format, as long as you give 77. Zeuss, D., Brandl, R., Brändle, M., Rahbek, C. & Brunzel, S. Global warming appropriate credit to the original author(s) and the source, provide a link to the Creative favours light-coloured insects in Europe. Nat. Commun. 5, 1–9 (2014). Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory Acknowledgements regulation or exceeds the permitted use, you will need to obtain permission directly from This work and Y.M. were supported by Marie Skłodowska Curie H2020-MSCA-IF-2017 the copyright holder. To view a copy of this license, visit http://creativecommons.org/ 795890 project EXTINCT of the European Commission; as well as, initially by Foun- licenses/by/4.0/. dation for Biodiversity Research and CESAB (Centre for the Synthesis and Analysis of Biodiversity, France) project LOLA. We thank the eBMS (https://butterfly- monitoring.net/) for the data and their volunteers for providing the data needed for the © The Author(s) 2022 study; and Roger Vila and the three reviewers for their revision on the manuscript. COMMUNICATIONS BIOLOGY | (2022)5:143 | https://doi.org/10.1038/s42003-022-03088-3 | www.nature.com/commsbio 9