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; Novella-Fernandez et al., 2023; 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 (Chen et al., 2011; Engelhardt et al., 2022; Grewe et al., 2013; 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).

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 result from 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 (Neate-Clegg et al., 2024; Stephens et al., 2016), including among insects (Kerr et al., 2015; Neate-Clegg et al., 2024; Pinkert et al., 2022), but potential relationships between phenological and geographic responses have not yet been investigated at continental scales. Functional 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 are inconsistent across taxa and regions, 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 animal species with high latitude ranges have been shown to exhibit smaller range shifts (MacLean and Beissinger, 2017; Pinkert et al., 2022), 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 a rare opportunity to investigate species’ spatiotemporal responses at larger taxonomic and spatial scales, particularly as most work has been done at national scales (Córdoba-Aguilar et al., 2023; Kalkman et al., 2018; Sandall et al., 2022). Due to their use of aquatic and terrestrial habitat across life different stages, dragonflies and damselflies are also considered indicator species for both terrestrial and aquatic insect responses to changing climates (Hassall, 2015; Pinkert et al., 2022; Šigutová et al., 2025), giving the study of these species broad relevance for conservation. 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; Rocha-Ortega et al., 2022; Suhonen et al., 2022), and rates of decline and expansion within limited geographic scopes (Powney et al., 2015; Rapacciuolo et al., 2017; Rocha-Ortega et al., 2020). 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, differences in habitat connectivity and dispersal ability may constrain range shifts for lentic species (those species that breed in slow moving water like lakes or ponds) and lotic species (those living in fast moving-water) in different ways (Kalkman et al., 2018). More southerly lentic species may expand their range boundaries more than lotic species, as species accustomed to ephemeral lentic habitats better dispersers (Grewe et al., 2013), yet lotic species have also been found to expand their ranges more often than lentic species, potentially due to the loss of lentic habitat in some areas (Bowler et al., 2021). While warm-adapted species with more equatorial distributions could expand their ranges poleward following warming (Devictor et al., 2008), they could also increase in abundance in this new range area relative to species that historically occupied those areas and are less heat-tolerant (Powney et al., 2015).

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 functional 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 species able to use both lentic and lotic habitats would shift their phenologies and geographies more than those able to use just one habitat type, 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). 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, 2010; Zattara and Aizen, 2021). 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. If a species was found on both continents, we only retained observations from the continent that was the most densely sampled. If we merged data for one species found on both continents, we could not perform a cross-continental comparison. However, if the same species on different continents was treated as different species, this would lead to uninterpretable outcomes (and the creation of pseudo-replication) in the context of phylogenetic analyses. In addition, species found on both continents did not have sufficient data to meet criteria for the phenology analysis.

We followed widely accepted methods to determine species range boundaries (Devictor et al., 2012, 2008; Kerr et al., 2015), although other methods exist that are appropriate for different data types and research questions i.e. (Ni and Vellend, 2021). We assigned species presences to 100×100 km quadrats, a scale that is large enough to maintain adequate sampling intensity but still relevant to conservation and policy (Soroye et al., 2020), 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 from Canada, the United States, and Northern Mexico, 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 1:

Figure 2:

Temperature Variability

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 main flight period for odonates in this study, for each year included in the historical and recent periods. We calculated the coefficient of variation 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; glm command in R) 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).

Range geographies and functional traits

To assemble trait data for the 76 species in the database, we used field guides (Cannings, 2002; Jones et al., 2008; Paulson, 2012) and existing trait databases (Powney et al., 2014; Waller et al., 2019). We considered any evolved morphological, physiological, behavioural, or life history characteristic as a functional trait (Beissinger and Riddell, 2021). Geographic range and associated climatic characteristics are often considered ecological traits, as they are consequences of functional traits and their interactions with geographic features (Bried and Rocha-Ortega, 2023; Chichorro et al., 2019). Such ecological variables may predict species’ responses to climate change, and can add significant value to predictive models (MacLean and Beissinger, 2017). We identified whether species’ ranges were more northern, southern, or both northern and southern (both), and determined the range size of each species by counting the number of quadrats occupied by that species in the historical time period.

Along with the geographic and climatic attributes, temperature variability and distribution, we selected four functional 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 (Cannings, 2002; Jones et al., 2008; Paulson 2012; Powney et al., 2014; Waller et al., 2019). 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. Egg laying habitat was assigned according to whether species use exophytic egg-laying habitat (i.e. eggs laid in water or on land, relatively larger in number), or endophytic egg-laying habitat (i.e. eggs laid inside plants, usually fewer in number); species using exophytic habitats are associated with greater northward range limit shifts (Angert et al., 2011). 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 the vast majority of 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 by modelling phenology shift as the dependent variable, and range shift and continent as independent variables. We used both species-level frequentist (GLM; glm function in R) and Bayesian (Markov Chain Monte Carlo generalized linear mixed model, MCMCglmm; Hadfield, 2010) models 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 molecular phylogenetic tree published by the Odonate Phenotypic Database (Waller et al., 2019), which used a morphological and taxonomic phylogeny as the backbone tree, allowing species to move within their named genera or families according to molecular evidence (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 functional 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. Non-significant variables, specifically all functional traits, were removed from the final geographic range shift model. No effects were significant in the model of phenology shifts. 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 (52 of 76) 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, although some species showed range retractions (Figure 4). Most species (41 of 66) maintained or advanced their emergence phenology (mean of −2.71 days in emergence phenology shifts among all species; Figure 4). 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.

The effect of species’ range shifts on phenology range shifts was significant in our model investigating the relationship between these responses, indicating that species shifting their northern range limits to higher latitudes also showed stronger advances in their emergence phenology (Figure 3). This result that was consistent in GLM and Bayesian analyses (p < 0.01; Table 1), and was maintained across North America and Europe, with no effect of continent in the model. This trend was consistent among both dragonflies and damselflies, although there was 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 (R² = 17% and 15%, for GLM and MCMCglmm models, respectively).

Figure 3:

Table 1:

Drivers of range and phenology shifts

Range geography and climate variability, rather than functional traits, predicted range shifts in both North America and Europe, with range geography being consistently the strongest predictor. Species’ functional traits did not relate to the extent of observed 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 species or species present in both the north and south (Table 2; p = 0.002 and 0.004, model R² = 26.6% and 23.7%, for GLM and MCMCglmm models, respectively), with no effect of range size on range shifts. Species experiencing smaller changes in interannual temperature variability also had a higher likelihood of northern range limit shifts (p = 0.0005 for GLM, p = 0.002 for MCMCglmm). Results from the GLM and MCMCglmm models were qualitatively similar, however, a smaller amount of model variance was explained when phylogeny was accounted for. Emergence phenology shifts were not affected by species’ traits, range geography, nor climate variability; due to this, model results are not displayed here. While range and phenology 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). Adding a phylogeny term to the MCMCglmm models also failed to produce a pattern different to the GLMs, and model performance did not improve when we accounted 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, but not functional 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 species or species present in both the north and south. 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 precursors to range expansions (Mair et al., 2014). Further mechanistic understanding of these processes requires abundance data. 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). Emerging mean conditions in areas adjacent to the ranges of southern species may offer opportunities for range expansions of these relative climate specialists, which can then tolerate climate warming in areas of range expansion better than more cool-adapted historical occupants (Day et al., 2018). Adaptive evolution and plasticity may enable higher population growth rates in newly-colonized areas (Angert et al., 2020; Usui et al., 2023), but this possibility can only be directly tested with long term population trend data. While some species experienced range retractions, these may result from sampling variability or stochastic population fluctuations along the northern range edge.

Increasing frequency and severity of extreme weather limited species’ geographical range responses (Table 2). This trend was independent of functional 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 these shifts in North American and European Odonata. 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 also helps to predict extinction risk during climate change (Kerr, 2020; Rocha-Ortega et al., 2020). The extent to which species’ traits actually determine rates of range and phenological shifts, rather than occasionally correlated with them, is worth considering further, but functional traits do not systematically drive patterns in these shifts among Odonates in North America and Europe.

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).

Given that species’ functional traits did not predict temporal or geographic responses, it is unsurprising that species’ responses were also independent of phylogenetic history (Franke et al., 2022). 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 species that co-occur may experience divergent selection within their climate tolerances (Schuetz et al., 2019).

It remains unclear if range and phenology shifts relate to trends in abundance, but our results suggest that there are clear ‘winners’ and ‘losers’ under climate change (Figure 4). 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). Species that did not shift their ranges northwards or advance their phenology 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). This group also includes Coenagrion resolutum, a common North American damselfly (Swaegers et al., 2014), for which we could not find evidence of decline. This may be due in part to the greater area of intact habitat available in North American compared to Europe, enabling C. resolutum to maintain larger populations that are less vulnerable to stochastic climate events. Still, this and other species failing to shift in range or phenology should be assessed for population health, as this species could be carrying an unobserved extinction debt. 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.

Figure 4:

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. 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. Climate change is associated with increasing variability as well as shifting mean conditions, contributing to species decline and even local extinction risks (Duffy et al., 2022). In this study, where species are found (i.e. their range geographies) determines whether they are exposed and respond to such negative pressures. Simultaneous consideration of shifts in range and phenology is a powerful and necessary approach to test aspects of species’ vulnerabilities to rapid global changes. By considering both the seasonal and range dynamics of species, emergent and convergent climate change responses across continents become clear for this well-studied group of predatory insects.

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).

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)

Additional information

Funding

Natural Sciences and Engineering Research Council

University of Ottawa