Species responses to weather anomalies depend on local adaptation and range position

Degree of local adaptation

We used count data from 34 butterfly species whose populations have been previously seen to show a clear response to specific climatic anomalies⁸ (Supplementary Figs. 1–34). These species have been previously classified as locally (N = 21) or globally (N = 13) adapted to the climatic regime based on their degree of local adaptation⁸. The latter is a continuum value measuring the species-specific sensitivity to climatic anomalies at two spatial scales: local (i.e., climatic deviations from the average across all years at that specific site) and global (i.e., climatic deviations from the average climate across all years and all sites where the species is found); quantified as the difference in the variance explained by the models:

where N_it is the annual abundance index calculated per site (i) and time (year, t); and W_it is the standardized climatic anomaly of the weather variable and for the temporal window (t’) most affecting the dynamics of the species at the site (i), either at the local or the global scale⁸. We included a density dependence term (logN_ijt-1) because it is important in butterfly populations^35,46, as well as the quadratic slope coefficient of the climatic anomaly to detect potential for non-linear responses. Site was set as a random intercept, except for species with a number of sites <10, and errors were set as normally distributed after verification.

The climatic variable W_it’ was calculated for temperature, precipitation and aridity (as a combination of temperature and precipitation) since butterflies, and many other organisms, can be affected by each of these variables^7,8,30,47. Likewise, the temporal window (t’) was defined for each life stage, given their relevance at modulating the effect of climate on population dynamics for many organisms, including butterflies⁴⁸. These were: pre-flight, flight, post-flight or overwintering period⁷ and year t and t-1 (to account for delayed responses), set per species and bioclimatic zone. Each period was defined using the annual flight curve distribution from the relative abundances per species and zone (i.e., using a generalised additive model fitted annual species phenology^34,49) as follows: flight period was set as the dates between the 10th and 90th percentiles of a species abundance distribution per year and bioclimatic zone for both univoltine and multivoltine species (i.e., species with more than one reproduction per year); pre-flight period was then set from February to the 10^th percentile of the species flight period; post-flight period from the 90th percentile to the end of October for year t-1 (year t was discarded because adults would be perished, so they have not possible effect on the counts or on the population dynamics of the species in that year); and overwintering period fixed from November to January (of year t – 1 to year t)^7,8.

The most sensitive scale of adaptation (local or global), weather variable and temporal window (phenological period and year t or t − 1; Supplementary Table 1) per species were then selected based on AIC best model fit in relation to population growth rate (Eq. 1 and Eq. 2). Comparison with null model accounting only for density dependence was previously tested⁸.

The degree of local adaptation was then quantified as the difference between the variance explained (R²) by the best model in relation to the anomalies at the local scale versus the variance explained by the model at the global scale⁸. Hence, values of the degree of local adaptation range [−1, 1] from most sensitive to anomalies at the entire species distributional scale (globally adapted species) to most sensitive to those anomalies at the local (site) scale (locally adapted species). The degree of local adaptation of the selected 34 species ranged from [0.004, 0.7] for local, [0.004, 0.06] if removing the outlier Parnassius apollo, with median 0.03 for both cases; and [−0.23, −0.002] for global species, [−0.07, −0.002] if removing the potential outliers Laeosopis roboris and Cupido osiris, median −0.02 for both cases (Supplementary Fig. 72 and Supplementary Table 1).

Data gathering

Species count data were collected via the long-term Butterfly Monitoring Schemes carried out in Finland, Spain and the UK covering six out of the ten bioclimatic regions across Europe (Supplementary Fig. 73). The schemes consist of a network of sites where volunteers perform weekly counts of butterflies along a set of transects following the standardized ‘Pollard Walk’ methodology (Pollard, 1977). Monitoring is done during the butterfly flight season, which varies depending on the climatic zone within the range of beginning of March to end of September. The three schemes differ in starting year and number of surveyed sites: Finland (1999, n_sites = 107), Spain (1994, n_sites = 130) and the U.K. (1976, n_sites = 2128). Therefore, we used Finland as the limiting country to set the range of study years (1999–2017). To have sufficient data per site, we only used data from those transects with at least ten years of interannual population change between 1999–2017, leading to 53 transects in Finland, 59 in Spain and 701 in UK. Counts were transformed to an index of abundance (N_jit) per species, site and year using generalized additive models that account for missing counts, and spatial and temporal variation in the species phenology^34,49. However, to assure consistent estimates, we excluded indices of abundance with more that 50% of weeks missing data. We used this annual species index of abundance to calculate the annual population change as log(N_jit/N_jit-1).

We extracted the site daily temperature (in degrees Celsius) and precipitation (in millimetres) at 0.1° spatial scale (~11 km) for each site from the European Climatic and Assessment Dataset project (ECAD)^50,51, to then calculate the climatic anomalies and to define the bioclimatic niche of the species (see Bioclimatic niche construction). Potential interactions between temperature and rainfall were accounted for by calculating the standardized aridity index^47,52:

where SAI stands for Standardized Aridity Index (hereafter aridity), T for mean temperature, P for total precipitation, sd for standard deviation calculated per site (i) and year (t).

For each species, we calculated the site climatic anomalies as:

where W_it’ is the anomaly, \(\bar{W}\)_it’ and \(\bar{W}\)_i are the mean temperature, the total precipitation or the aridity index (SAI) over the time series, depending on the climatic variable to which each species was most sensitive to, calculated per site (i) and, for \(\bar{W}\)_it’, per temporal window (t’; phenological period and year t or t-1) during which the climatic anomalies most affected the population dynamics of the species (Supplementary Table 1)⁸.

Bioclimatic niche construction

We defined the species-specific bioclimatic niche following the realised niche concept, which delimitates the range in which a species can persist based on the physical environmental conditions in combination with biotic constraints (Hutchinson, 1957). For each species, the bioclimatic niche was constructed based on the species distribution along the minimum-to-maximum climatic values where the species was detected during the study period along all studied sites; using, per species, the climatic variable most affecting its population responses (temperature, precipitation or aridity; Supplementary Table 1).

Each species bioclimatic niche was calculated taking the annual mean climatic value of the temporal window associated to the species, at the site where the species was recorded as:

Where W_it’ is the mean climatic value per site (i) and temporal window (t’) and W is the observations of climatic values from all sites and years¹⁵. The bioclimatic niche for precipitation was calculated with its inverse values, so that the bioclimatic niche was coincident with the species distribution along the temperature and the aridity index (e.g., coldest marginal edge will also have higher levels of precipitation).

This bioclimatic niche calculation preserves the relative difference between the climatic values of the sites where a species was observed but standardises the climatic range differences between species. As a result, the minimum mean climatic value each species was ever detected was given a score of −1 and the maximum of 1. Thus, every site and species combination had a position within a range between −1 and 1, delineating the bioclimate leading and trailing margins, respectively and broadly corresponding with the edges of the species distribution (Supplementary Fig. 35).

To assure the alignment of the bioclimatic niche and the species population range position, we performed species specific Pearson correlations between them (Supplementary Fig. 74). To do so we used the relative range position, defined as the relative position of each population within the latitudinal distribution of its species and calculated based on the formula⁷:

where L_i is the latitude, being \(\bar{{L}_{i}}\) the average latitude of the species latitudinal distribution set as a vector, weighted by the number of data points (species abundance) from each site. RPP is then expressed as a proportional range position standardised from 0–1 for each species, with 0 zero relating for lower latitudinal position (hence trailing marginal edge of the species) and one for higher latitudinal position (hence leading marginal edge of the species). While we lack of the entire distributional information for each species, this method allowed us to test for a proxy of correlations in latitudinal gradients, hence accounting for trailing and marginal edge margins.

Statistics and reproducibility

We fitted linear mixed models to test if the variation of the species populations responses, in terms of population change, to the climatic anomalies along their bioclimatic niche differed between globally and locally adapted species. We set population change (log(N_jit/N_jit-1)) as the response variable, and fitted local climatic anomaly (at the site level) in an interaction with the position of the site within the species bioclimatic niche. We also included density dependence (logN_ijt-1) and the quadratic slope coefficient of the climatic anomaly to account for possible non-linear responses for each species. Species and site were set as random effects, and the error as normally distributed:

where N_jit is the annual abundance index of the species j at the site i and year t, W_it’ is the climatic anomaly variable at the site i and temporal window of the species (phenological period and year t or t-1; Supplementary Table 1) at the local scale, and B_ji is the bioclimatic position of the population of species j at the site i. Models were fitted with the species most affecting standardised climatic variable W (temperature, precipitation and aridity), and temporal window⁸. Climatic anomalies were standardised separately per variable (temperature, precipitation and aridity), and precipitation was used in inverse (i.e., multiplied by −1). Models were fitted pooling all species together but separately for locally (degree of local adaptation] 0, 1]) and globally (degree of local adaptation] 0, −1]) adapted categories. As a post-hoc analysis we also fitted the models with a phylogenetic linear mixed model for locally and globally adapted species separately, to account for phylogenetic non-independence across species because we previously found a phylogenetic signal in the degree of local adaptation⁸ (Supplementary Table 4).

We also fitted linear mixed models to test whether population trends (i.e., abundances over time) were more stable for locally than for globally adapted species, and test the relation between the trends and the species bioclimatic niche. Models were fitted with population abundance (set as the continuous index of abundance N_jit) as the response variable, transformed as log(N_jit + (N_jit² + 1)^1/2) to account for zero-values, and the cumulative numbers of years of observations per species in interaction with the population position within the species bioclimatic niche as explanatory variables. We used the transformation log(N_jit + (N_jit² + 1)^1/2) rather than direct logarithmic as it performs similarly to a logarithmic transformation but zero remains defined avoiding the need to remove these data or add an arbitrary value⁵³. The cumulative numbers of years set per species (nested random effects) and site were set as random effects and the error as normally distributed:

where N_ji is the annual abundance index of the species j at the site i and time t, Z_jt is the cumulative number of years of observations of the species j at the time t, and B_ji is the bioclimatic position of the population of species j at the site i.

To increase robustness of all our analyses we also performed separate models with alternative stricter definitions for local and global adaptation of species and without potential outliers (Supplementary Fig. 71). For locally adapted species, we modelled population change and abundance trends for those species with degree of local adaptation constrained to [0.025, 1] and without Parnassius apollo, and for globally adapted species from [−1, −0.025] and without Laeosopis roboris and Cupido Osiris. This led to four models for each species category (local and global; Supplementary Table 2 and Supplementary Table 5). We did not model population change using more conservative values of the degree of local adaptation because a more conservative approach reduced the number of species to <10.

We did a stepwise model reduction for interaction or additive effect between the climatic anomalies and the population position within the species bioclimatic niche. We also did model reduction between quadratic and linear response, using the full model to account for potential broadly linear responses of globally adapted species (all AIC model values provided in Supplementary Tables 6 and 7). We conducted all analyses in R 3.6.1 (R Core Team, 2019), using the package lme4⁵⁴ to fit our models and the package MuMin for model averaging⁵⁵.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.