Introduction

Climate change alters climatic means and increases the frequency of extreme weather events, exposing species to conditions outside of their tolerances, and often leading to population declines (Goulson, 2019; IPCC, 2021). Species may avoid extreme conditions by dispersing to new areas where conditions pose fewer weather-related challenges, often leading to poleward range expansion (Davis et al., 2005; Lawlor et al., 2024). Species’ biological timing could also shift, through adaptation or phenotypic plasticity, with earlier warming advancing the timing of early-season activities and life-history events (Davis et al., 2005; Hällfors et al., 2024; Parmesan and Yohe, 2003). Both species’ geographical ranges and seasonal timing depend strongly on climate and habitat conditions, with shifts in space and time permitting species to remain within the limits of their ecological niches (Figure 1; Chen et al., 2011; Menzel et al., 2006; Parmesan and Yohe, 2003). This allows populations to grow, despite changing environments, and reduces the risk of climate debt and extinction (Devictor et al., 2012; Franks et al., 2018; Lustenhouwer et al., 2018; Saino et al., 2010; Souza et al., 2023; Urban, 2015).

Figure 1:

Positive population trends can be stronger in species that shift both their range and phenology (Hällfors et al., 2024). Greater phenological plasticity under warmer spring temperatures may increase reproductive success, leading to greater population growth and range expansions (Macgregor et al., 2019), as positive or stable trends in species abundance and habitat availability are essential for range shifts (Mair et al., 2014; Platts et al., 2019). However, there could also be a trade-off between phenological and geographical shifts (Amano et al., 2014; Hassall, 2015; Socolar et al., 2017). Species with greater dispersal abilities may have less need for phenological shifts as they track their climatic niche through space, while weaker dispersers may be confronted with greater selective pressure to shift phenology within their range (Amano et al., 2014; Hassall, 2015; Socolar et al., 2017). Since range shifts can also be driven by extirpations at species’ trailing range edge (Parmesan et al., 1999), greater phenological shifts may mitigate the need for range shifts, as species better tolerate new climatic conditions. Cross continental studies report converging effects of climate change on species’ range shifts and abundances (Kerr et al., 2015; Stephens et al., 2016), but potential relationships between phenological and geographic responses have not yet been investigated at continental scales.

Traits, such as dispersal, may determine species’ spatial and temporal responses to climate change (Chen et al., 2011; Kharouba et al., 2009; Schuetz et al., 2019; Zografou et al., 2021). However, these relationships generalize poorly taxonomically and geographically, and cross-continental tests have not been attempted (Angert et al., 2011; Buckley and Kingsolver, 2012; Estrada et al., 2016; MacLean and Beissinger, 2017).

Geographic locations and environmental characteristics of species’ ranges may also predict range shifts, as species with high latitude ranges have been shown to exhibit smaller range shifts (MacLean and Beissinger, 2017), while increasing local temperature and loss of natural landcover may drive range retractions (Pacifici et al., 2020). However, exposure to extreme climate events (e.g. drought, heat waves, or storms) within species ranges may disrupt species’ dispersal abilities and capacities to tolerate new conditions (Kerr, 2020; Román-Palacios and Wiens, 2020). Exposure to thermal anomalies can rapidly change entire communities and create shifts towards new ecosystems, sometimes leading to local declines (Day et al., 2018; Grant et al., 2017; Harris et al., 2018; Román-Palacios and Wiens, 2020).

While global change research on insects often emphasizes butterfly and bee taxa, recently assembled databases of odonate observations provide opportunities to answer new questions (Kalkman et al., 2018). Dragonflies and damselflies respond rapidly to climate change, both spatially (Ball-Damerow et al., 2015; Grewe et al., 2013; Hickling et al., 2005; Powney et al., 2015; Rapacciuolo et al., 2017; Sirois-Delisle and Kerr, 2022) and temporally (Hassall, 2015; Hassall et al., 2007; McCauley et al., 2018), making them an excellent indicator species (Hassall, 2015). There is some evidence that functional traits relate to odonates’ interspecific variation in range shifts (Angert et al., 2011; Grewe et al., 2013), phenology shifts (Diamond et al., 2011; Gutiérrez and Wilson, 2021; Zografou et al., 2021), extinction risks (Cardillo et al., 2008; Cooper et al., 2008; Fritz et al., 2009; Suhonen et al., 2022), and rates of decline and expansion within limited geographic scopes (Powney et al., 2015; Rapacciuolo et al., 2017). While relationships between morphological traits and range boundaries have been shown for some groups (i.e. Rundle et al., 2007), these may depend on species’ geographic context. For example, lentic species (i.e. species breeding in slow moving water like lakes or ponds) tend to have larger ranges and higher dispersal capacities than lotic species (i.e. fast moving water-dwelling species), and may track their climate niches more closely (Grewe et al., 2013; Hof et al., 2006). Lotic habitat may be more difficult to reach, and lotic species may face stronger anthropogenic pressures in places where their diversity is high (Kalkman et al., 2018), affecting species abilities to track changes in climatic seasonal timing. Lentic species have been shown to expand their range boundaries more than lotic species, but only among species with more southern distributions (Grewe et al., 2013). Abundances of European odonates with ranges extending south from sample landscapes have been found to respond more positively to climate change than species with ranges extending north from sample locations, due to their geographic position and thermal tolerance (Powney et al., 2015). Species previously limited by tolerance to low temperatures could expand their ranges following warming (Devictor et al., 2008).

In this study, we tested whether species with stronger geographical range shifts also advanced their emergence phenologies, or if one response offsets the need for the other. We also asked whether biological traits, range geography (i.e. southerly vs. northerly), or temperature variability predict range shifts at species’ northern range limits, and whether these factors can also predict shifts in species emergence phenology. We predicted that species would exhibit shifts in both geography and phenology, as we expect species shifting in phenology to have larger populations, increasing likelihood of dispersal and successful range shifts. We predicted that generalist species would shift their phenologies and geographies, as generalists outperform specialists as climate and land uses change (Ball-Damerow et al., 2015, 2014; Hassall and Thompson, 2008; Powney et al., 2015; Rapacciuolo et al., 2017); in contrast, we expect specialists to fail to shift in both space and time, making them more vulnerable to extirpation and extinction. Alternatively, species might respond to rapid climate changes in ways that reflect their geographical position, indicating that where they are found is a better predictor of their conservation risk from global change than their intrinsic biological characteristics.

Methods

Biological records

We assembled ~2 million observations records for North American and European odonate species collected between 1980 and 2018. Data sources included online dataset aggregators GBIF (http://gbif.org/) and Canadensys (http://www.canadensys.net/), Odonata Central (Abbott, 2020), and other institutions (see Acknowledgments). While odonates were sampled opportunistically, biases associated with data that are not systematically collected are less likely to affect trends at large spatial and temporal scales, particularly if data are obtained from multiple independent sources (Pyke and Ehrlich, 2010a; Zattara and Aizen, 2021a). We removed records with incomplete or missing species identification, year, or locality information. We selected unique observations for species, location, year, and Julian day of collection, and restricted the data to continental North America and Europe. We mapped species-specific observations using ArcGIS software (ESRI, 2019), and qualitatively verified species ranges. In the event that the same species was found in both Europe and North America, we selected records from the continent containing the largest number of observations (Supplementary Information).

We assigned species presences to 100 × 100 km quadrats to identify the best sampled species. We excluded species found in fewer than 50 quadrats to increase the likelihood of accurately predicting the position of species’ northern range boundaries. We retained ~1,100,000 records comprising 76 species (Figure 2). Observation records were separated into two time periods, to compare species’ recent phenologies and northern range limits (2008 to 2018) to conditions in a historical time period (1980 to 2002).

Figure 2:

Climate data

Temperature variability during a species’ flight season may impact its ability to establish in new locations, or shift its emergence timing the following year. We downloaded a high-resolution gridded dataset for monthly average daily maximum temperature from the Climatic Research Unit (New et al., 2002). We extracted average values per 100 × 100 km quadrat across the months of April to October, covering the flight period for odonates in this study, for each year included in the historical and recent periods. We calculated the standard deviation per quadrat for each time period and averaged these values per species and time period to measure interannual temperature variability during species’ flight season across their range. We used the difference between recent and historical measures as an estimate of change in temperature variability.

Spatial and temporal change metrics

To limit potential effects of temporal and spatial biases, we generated range and phenology shift metrics with specific criteria for quadrats and species selection (Bartomeus et al., 2018; Gaiji et al., 2013; García-Roselló et al., 2015); we retained 76 species for range shift estimates (Supplementary Information). Species’ northern range boundaries were calculated using the mean of the 10 most northern points of each species range in both the historic (1980 – 2002) and recent (2008 – 2018) time periods, measured in kilometers from the equator (as in Kerr et al., 2015). We used the difference between range limit positions in the historic recent time periods to estimate species’ northern range limit shifts.

Since spatiotemporal biases sometimes inflate range shift measurements (Kujala et al., 2013), we used null models to test whether observed range shift estimates were robust, and the extent to which those responses differed from expectations arising because of rising sampling intensity over time. 1000 randomized datasets were created with the same number of species-specific geographical points per time period as the number of actual observations. Maximum and minimum latitude and longitude values were held constant relative to observed values. Northern range limit shifts were calculated as the difference in the mean of the 10 most northern points between the historic and recent time periods. We used generalized linear models (GLMs) to test whether species-specific range limit shifts in each iteration predicted range shifts measured from the observation data (Figure S1, S2). We found that observed northern range limit shifts are not consistent with expectations derived from changes in sampling intensity.

We estimated species-specific emergence phenology for each time period in 200 × 200 km quadrats; using larger quadrats increases probabilities of detecting signals of emergence phenology, which may otherwise be lost due to gaps in data density. We retained quadrats that contained at least 25 observations for a given species in both time periods. To estimate phenology per area, we used the weib.limit function of the WeibullR R package (Pearse et al., 2017). This function uses the Weibull distribution to estimate the Julian day of a species’ first appearance, and is especially useful to measure the timing of phenological events in sparsely sampled datasets. Techniques such as using the average of the n^th first observations of phenological events, or the n^th percentile flight dates (Brooks et al., 2014; Robbirt et al., 2011), tend to overestimate the timing of biological events due to temporal bias towards later days in the species’ active period (Pearse et al., 2017). We retained emergence estimates between March 1^st and September 1^st, as well as species and quadrats that showed a difference in emergence phenology of −25 to 25 days, −30 to 30 days, or −35 to 35 days between both time periods, to include only phenology shifts that could be biologically meaningful to environmental climate change. 68 species found across 63 quadrats met these criteria. Large changes in phenology are likely explained by other anthropogenic or natural factors, or could occur due to noise in the data, since these phenology calculations per region are extremely data intensive. We calculated the difference in the day of emergence per quadrat between both time periods, as well as mean phenology change across all quadrats for each species. The number of quadrats per species used to calculate their mean phenological shift varied between 2 and 46.

We used null models to assess whether our approach to estimating phenology shifts was robust, as potential issues may arise due to spatiotemporal biases in the underlying data (Kujala et al., 2013). We constructed 1000 randomized datasets of species’ hypothetical days of occurrence, using the same number of quadrats, and observations within quadrats, as in each time period of the observation records. We assigned the maximum and minimum Julian day of occurrence from the observation records to limit values in the randomized datasets. We applied the same method and criteria of inclusion to the randomized datasets as we did to measure phenology shifts from the observation data. GLMs were built to test whether phenology shifts calculated using 1000 random datasets predicted the phenology shifts that we measured. No discernable pattern emerged, indicating that observed shifts in phenology are not consistent with expectations derived from differences in sampling intensity over time (Figure S3, S4).

Ecological traits

To assemble trait data for the 76 species in the database, we used field guides (R.A. Cannings, 2002; C.D. Jones et al., 2008; D. Paulson, 2012) and existing trait databases (G. D. Powney et al., 2014; Waller et al., 2019a). We considered any evolved morphological, physiological, behavioural, or life history characteristic as a trait (Beissinger and Riddell, 2021). Geographic or climatic characteristics (e.g. range size, geographic position of species’ ranges, climatic conditions of the range) were not considered traits, as they are not products of evolution, but rather consequences of environmental, ecological, biogeographical, or evolutionary patterns linked to range shifting mechanisms (Beissinger and Riddell, 2021). Such variables may predict species’ responses to climate change nonetheless, and can add significant value to predictive models (MacLean and Beissinger, 2017).

Along with two geographic and climatic attributes, including temperature variability and distribution region (northern, southern, or widespread), we selected four traits likely to be biologically relevant to spatial and temporal responses to climate change: flight duration, breeding habitat type (lotic, lentic, or both), egg laying habitat (exophytic vs. endophytic), and body size. Species’ flight period was measured as the total number of days of the flight period, estimated from the approximate time of the month of average first and last appearances. Breeding habitat was assigned according to a species’ uses of lotic, lentic, or both habitat types. Body size corresponded to the mean length of the abdomen of each species. We excluded overwintering stage and range size from our analysis as data were incomplete for many species, and excluded migratory behavior as nearly all species included in the study were non-migratory.

We tested for correlations amongst all predictors by calculating the Predictive Power Statistic (PPS) and Pearson correlations among traits (van der Laken, 2021): we found no evidence of correlation.

Statistical analyses

We conducted statistical analyses using R Statistical Software (R Core Team, 2019). All continuous variables were transformed into Z-scores using the scale function in R. First, we investigated whether there was a relationship between species’ range and phenological shifts. We used a by-species GLM and a Markov chain Monte Carlo generalized linear mixed model (MCMCglmm; Hadfield, 2010); we compared frequentist and Bayesian approaches to improve the robustness of the results. We included a term to account for phylogeny in the MCMCglmm model, as species that are closely related are likely to have similar traits. We used the phylogenetic tree published by the Odonate Phenotypic Database (Waller et al., 2019a), which was generated using a combination of DNA-sequences and species’ taxonomy based on morphology (Waller and Svensson, 2017). Trace and density plots for the MCMCglmm model revealed no issues with autocorrelation or model convergence (Figure S7).

Next, we investigated whether biological traits, range geography, or temperature variability predicted range shifts at species’ northern range limits, and whether the same predictors explaining range expansions could also predict shifts in species emergence phenology. We constructed two sets of GLMs, in addition to two sets of MCMCglmms accounting for phylogeny; one of each with changes in species’ northern range limits as the response variable, and the other with changes in emergence phenology as the response variable. We included species traits, species range geography, and temperature variability experienced by species as predictor variables. Trace and density plots for the MCMCglmm models did not indicate limitations related to autocorrelation or model convergence (Figure S8).

In addition to the inclusion of phylogeny in statistical models to account for potential data non-independence, we measured the phylogenetic signal in range and phenological shifts. We used the phylosig function of the Phytools package version 0.7-70 (Revell, 2012), which calculates phylogenetic signal using Pagel’s lambda and Blomberg’s K.

Results

Relationship between range and phenology shifts

Most species (69% of those included) expanded their northern range limits toward higher latitudes (mean range expansions of 180 km). The average range expansion across all species was 63 km northward (Figure S5). Some species showed range retractions (i.e. Libeullula luctuosa); small range retractions may result from sampling variability or stochastic population fluctuations along the northern range edge. Most species (60% of assessed species) maintained or advanced their emergence phenology (mean of −2.71 days in emergence phenology shifts among all species; Figure S6). Fewer species were included in phenology analyses than analyses of range shifts due to the data intensity required to capture phenology shifts at a study site (N=66 vs. N = 76). Many species (50%) showed both advancing emergence phenology and range expansions, while 10% of species showed neither range nor phenological shifts relative to historical baselines. Species that did not respond to climate change temporally or spatially included Coenagrion mercuriale, a European species that is listed as near threatened by the IUCN Red List (IUCN, 2021), and is projected to lose 68% of its range by 2035 (Jaeschke et al., 2013).

Species that shifted their northern range limits to higher latitudes also showed stronger advances in their emergence phenology (Figure 3), a result that was consistent in GLM and Bayesian analyses (p < 0.05; Table 1; R² = 17% and 15%, for GLM and MCMCglmm models, respectively). This trend converged across continents, and is consistent among both dragonflies and damselflies, although there is considerable interspecific variation in the magnitude of spatial and temporal shifts. Accounting for phylogeny did not improve model predictions and did not explain greater model variance (Table 1).

Figure 3:

Table 1:

Drivers of range and phenology shifts

Range geography and climate variability, rather than ecological traits, predicted range shifts in both North America and Europe, with range geography consistently the strongest predictor. Species’ ecological traits did not relate to the extent of observed phenological or geographical range shifts in tests using GLM and MCMCglmm models. Species with more southern distributions shifted their northern range limits towards higher latitudes more than northern and widespread species (p < 0.05; Table 2; R² = 23% and 20%, for GLM and MCMCglmm models, respectively). Species experiencing smaller changes in interannual temperature variability also had a higher likelihood of northern range limit shifts (p < 0.05). Results from the GLM and MCMCglmm models were qualitatively similar, however a smaller amount of model variance was explained when phylogeny was accounted for. Neither range geography nor climate variability predicted emergence phenology shifts. Even though both response types are related, this suggests that the mechanisms underlying phenological shifts are different than those underlying range shifts.

Table 2:

A phylogenetic signal may indicate that there are traits that determine species’ spatial and temporal responses to changing climate that were not measured in this study. Yet we detected no phylogenetic signal using Pagel’s lambda or Blomberg’s K in either geographical range or phenological responses (Table 2), despite pronounced recent warming. Adding a phylogeny term to the MCMCglmm models also failed to produce a pattern different to the GLMs, and model performance was not improved when accounting for phylogeny in our assessment of northern range shifts.

Discussion

In one of the first studies to investigate both shifts of phenology and range at a continental scale, we find that dragonfly and damselfly species show pronounced geographical and phenological shifts that converged across Europe and North America. Species expanding their ranges poleward also emerged earlier in the spring on both continents (Figure 3), with shifts predicted by range geography and climate variability, not ecological traits. These results suggest that some species are ‘universal winners’ with respect to climate change: they demonstrate the flexibility to respond both temporally and spatially to the onset of rapid climate change. Conversely, species that show neither geographic nor phenological shifts may be particularly vulnerable to climate change.

We found no evidence for a tradeoff between range and phenology shifts; instead, half of species shifted both range and phenology. Earlier seasonal timing allows species to stay within their climatic limits. This may lead to higher local population growth (Macgregor et al., 2019), although earlier emergence could expose individuals to early season weather extremes (McCauley et al., 2018). As only a small proportion of odonate adults undertake long range dispersal (Conrad et al., 1999), greater local population sizes should contribute to higher dispersal rates (Mair et al., 2014), facilitating range shifts (Kerr, 2020; Leroux et al., 2013). This is consistent with results from other taxa: among British butterflies, early emergence increased population growth and facilitated range shifts for species with multiple generations per year (Macgregor et al., 2019); Finnish butterfly species with the greatest population growth rates shifted both their phenology and ranges (Hällfors et al., 2021). Such population growth or maintenance, and therefore the potential for range shifts, is only possible if habitat is available (Mair et al., 2014). Future work should consider habitat availability alongside range and phenology shifts, as it may help explain why some species are able to shift their phenology but not their range.

Southern species were more likely to expand their ranges northward than northern or widespread species. Species’ ability to maintain large populations may be impaired in northern latitudes, where rates of climate change are high (IPCC, 2021), hindering dispersal and colonization that are necessary for range expansions (Mair et al., 2014). Southern species may have narrower niche breadths than widespread or northern species, and may respond more rapidly to climate change to track this narrower niche (Hällfors et al., 2024). As more southern species expand northwards they may be better able to tolerate local temperature extremes due to their warmer thermal niches, and should be better able to maintain populations than species with northern distributions (Day et al., 2018), facilitating further range shifts.

Increasing frequency and severity of extreme weather limited species’ geographical range responses (Table 2). This trend was independent of traits that are mechanistically linked to species’ climate change responses, such as dispersal ability or habitat preference. Extreme temperatures can reduce population sizes, leading to local extinctions (Román-Palacios and Wiens, 2020), and reducing the likelihood of range expansions (Mair et al., 2014). In odonates, experimental evidence has demonstrated that larval mortality rises with short-term extreme weather (Mccauley et al., 2015). Individuals that shift phenologies earlier in the season to avoid climate extremes could still be exposed to harmful conditions (Iler et al., 2021); for example, odonate populations that respond to unusually warm spring temperatures may experience high mortality if temperatures return to seasonal conditions. Species that experience extreme conditions may then be unable to successfully shift in time, reducing population sizes, and reducing the likelihood of range shifts.

In contrast to previous work demonstrating that range and phenology shifts are at least partially determined by species traits (i.e. Sunday et al., 2015; Zografou et al., 2021), no functional trait, or combination of traits, explained species’ shifts. While we could not capture all functional traits in this analysis, our results are consistent with other work that identifies climate velocity and sensitivity as the best predictors of range shifts and thermal preferences tracking in marine systems (Pinsky et al., 2013; Schuetz et al., 2019). Species’ tolerances to increasingly variable temperatures, rather than traits, also help predict extinction risk during climate change (Kerr, 2020). More work is needed to determine if traits consistently fail to predict range shifts, particularly at large spatial scales, although reliable observations at these scales are often lacking. Detailed mechanistic models that incorporate multiple independent processes may help untangle these patterns, although challenges with parameterization remain (Urban et al., 2022).

The geographic positions of species’ ranges determine the local pressures and environmental factors to which they are exposed (MacLean and Beissinger, 2017; Pacifici et al., 2020), potentially masking the effects of traits (Schuetz et al., 2019). This process could cause trait-related trends to differ across levels of biological organization (Srivastava et al., 2021), from local populations (where traits might be critical) to biogeographical extents (where traits might be unrelated to range or phenological shifts; Grewe et al., 2013; Gutiérrez and Wilson, 2021; Sunday et al., 2015; Zografou et al., 2021). Tests of local adaptation in critical traits, like thermal tolerances among larval life stages, would be valuable.

Given that species’ traits did not predict temporal or geographic responses, it is unsurprising that species’ responses were also independent of phylogenetic history. The phylogenetic approach did not improve model predictions in any model that we tested, and there was no phylogenetic signal in either response according to Pagel’s lambda and Blomberg’s K (Table 2). These results are consistent with previous work that found no phylogenetic trend in local odonate population extinctions (Suhonen et al., 2022). There may be strong variation in thermal niches among closely related species: species that are geographically isolated adapt to different local climates, while populations that coexist may experience divergent selection within their climate tolerances (Schuetz et al., 2019).

While more work is needed to determine if the range and phenology shifts reported here are associated with increases in species’ abundances, our results suggest that there are clear ‘winners’ and ‘losers’ under climate change. Climate ‘winners’, species that are shifting in space and time, may require more limited conservation intervention. Species expanding their ranges could be better supported if habitat area and connectivity are conserved, facilitating climate-driven range shifts (Littlefield et al., 2019). Species only shifting their phenologies may require further study, as phenology shifts may have positive or negative impacts on abundance (Iler et al., 2021). Climate ‘losers’, species that are failing to shift in both space and time, may require more direct conservation intervention, such as managed relocation (Richardson et al., 2009). Our analysis of phenology and range shifts should be repeated in other taxa, as it may offer a method of identifying conservation actions among species groups.

Here, we showed that odonate species exhibit convergent responses of range and phenology shifts across continents. While species with southern distributions were more likely to shift their ranges, increasing temperature variation limited geographical range responses among species in both Europe and North America. Simultaneous consideration of shifts in range and phenology is a powerful and necessary approach to measure species’ apparent vulnerability to rapid global changes. Understanding how range and phenology shifts vary across species, and what drives this variation, is increasingly urgent as climate change alters local and regional environmental conditions. This study gives us new insight into patterns of range and phenology shifts, as well as suggesting processes behind the patterns, offering a valuable method of understanding species’ vulnerability to change that should be explored in other taxa.

Supplementary Information

Odonate data

The full raw dataset includes ~2 million Odonata occurrence records. We first removed inadequate records from our primary dataset. Inadequate records were identified as those with incomplete information for species identification, year, or locality, inaccurate georeferenced points, and duplicate records. We retained ~1,100,000 unique species-location-date observations. There are 452,929 records in the first time period (1980 – 2002) and 641,590 observations in the second time period (2008 – 2018).

If a species was found on both continents, we only kept observations from the continent that was the most densely sampled. This simplified analyses in terms of calculating range and phenology shifts. It would not be logical to perform a cross-continental comparison when data points for a single species exist in both places and are merged together to create a single measurement of range shifts or location-specific measures of phenological change. This problem cannot be solved by treating species as different units because that leads to uninterpretable outcomes (and creation of pseudo-replication) in the context of phylogenetic analyses. Most importantly, it was not possible to calculate phenology shifts for both continents in these species because of the data intensity needed to make the calculations and the strict criteria we set for inclusion of species and quadrats in the analysis. Future work may focus on species found on both continents if the data and statistical analyses permit.

Phenology and range position estimates

We have not tried to interpret noise in phenology shift estimates, which could occur due to sampling issues, or to another biologically meaningful eco-evolutionary response. Our methods do not enable tests of those ideas, however. For example, Aeshna umbrosa has a surprisingly high phenology shift (>25 days) but a very strong range retraction (<-300 km). For this species, we retained 8 quadrats to calculate mean phenology shift estimate. It is possible that a section of the range was lost in which emergence was especially early, compared to other regions, due to local adaptations. This pattern may affect results where species appear to shift emergence dates earlier/later, but phenology shifts are affected by parts of the range that remain occupied.

We put in place criteria to make sure to include species with range shifts likely to result from climate change effects, rather than from sampling issues or to other anthropogenic pressures such as land use change, or sudden land use intensification. Nehalennia irene was thus removed, as its range positions were highly unusual, and likely due to other factors than global change (>800 km). Libellula quadrimaculata was removed due to sampling intensity discrepancies between time periods, having over 10,000 extra points in the second period compared to the first.

Temporal and spatial bias are likely to be present in opportunistic data, but they are less likely to impact long term factors of species’ distributions if including data from as many sources as possible, and that span across large geographic and temporal scopes (Pyke and Ehrlich, 2010; Zattara and Aizen, 2021). Here, we put in place several criteria in careful interpretation of results, as described in the Methods section of the main text. Further, among preliminary examinations of the data, we test for the effect of sampling onto our phenology and range shift estimates. We found no effect of sampling on range of phenology responses that we calculated here. We tested a linear mixed model with species as a random effect that assesses whether sampling intensity change per quadrat (where sampling intensity corresponds to sampling per species/quadrat) predicts species’ phenology change per quadrat. The p-value of the predictor variable, sampling intensity change, was not significant. Secondly, we built a Generalized Linear Model to assess whether greater sampling overall (measured as the number of quadrats in which each species is sampled) predicts species’ range shifts that we computed. The p-value of the predictor variable, sampling change per species, was not significant. If spatial bias was present here, we may expect the number of quadrats in which species are found to increase the position of species’ northern range boundaries.

Model information and statements

We report the full model statements below as executed with the MCMCglmm R package, including the priors, and the model of trait evolution to account for phylogeny. Prior structure followed standard practice according to continuous or binary response variable. The family “gaussian” or “threshold” was selected depending on the continuous or binary nature of the data, respectively. The number of iterations, burnin, and thinning parameters were tested and confirmed with a visual assessment of model convergence using the trace and density plots. It is common practice for the burnin parameter to correspond to 10% of the number of iterations.

Model 1

Aphylo ← vcv(phylogeny, model = “Brownian”, corr = T)

Ainv ← inverseA(phylogeny, nodes = “TIPS”, scale = F)$Ainv

prior ← list(R=list(V=1,nu=0.002))

MCMCglmm(Phenological shifts ~ Range shifts + Continent, ginverse = list(species = Ainv), family = “gaussian”, prior = prior, DIC = T, nitt = 150,000, thin = 100, burnin = 15,000, data)

Model 2

Aphylo ← vcv(phylogeny, model = “Brownian”, corr = T)

Ainv ← inverseA(phylogeny, nodes = “TIPS”, scale = F)$Ainv

prior ← list(R=list(V=1,nu=0.002, fix = 1))

MCMCglmm(Range shifts ~ Breeding habitat type + Distribution region + Flight duration + Egg laying habitat + Body size + Temperature variability, ginverse = list(species = Ainv), family = “threshold”, prior = prior, DIC = T, nitt = 150,000, thin = 100, burnin = 15,000, data)

Model 3

Aphylo ← vcv(phylogeny, model = “Brownian”, corr = T)

Ainv ← inverseA(phylogeny, nodes = “TIPS”, scale = F)$Ainv

prior ← list(R=list(V=1,nu=0.002))

MCMCglmm(Phenology shift ~ Breeding habitat type + Distribution region + Flight duration + Egg laying habitat + Body size + Temperature variability, ginverse = list(species = Ainv), family = “gaussian”, prior = prior, DIC = T, nitt = 150,000, thin = 100, burnin = 15,000, data)

Table S1:

Table S2:

Table S3:

Table S4:

Figure S1:

Figure S2:

Figure S3:

Figure S4:

Figure S5:

Figure S6:

Figure S7:

Figure S8: