Species are often expected to shift their distributions either poleward or upslope to evade warming climates and colonize new suitable climatic niches. However, from 18 years of fixed transect monitoring data on 88 species of butterfly in the midwestern United States, we show that butterflies are shifting their centroids in all directions, except towards the region that is warming the fastest (southeast). Butterflies shifted their centroids at a mean rate of 4.87 km yr-1. The rate of centroid shift was significantly associated with local climate change velocity (temperature by precipitation interaction), but not with mean climate change velocity throughout the species’ ranges. Species tended to shift their centroids at a faster rate towards regions that are warming at slower velocities but increasing in precipitation velocity. Surprisingly, species’ thermal niche breadth (range of climates butterflies experience throughout their distribution) and wingspan (often used as a metric for dispersal capability) were not correlated with the rate at which species shifted their ranges. We observed a high phylogenetic signal in the direction species shifted their centroids. However, we found no phylogenetic signal in the rate species shifted their centroids, suggesting less conserved processes determine the rate of range shift than the direction species shift their ranges. This research shows important signatures of multidirectional range shifts (latitudinal and longitudinal) and uniquely shows that local climate change velocities are more important in driving range shifts than the mean climate change velocity throughout a species’ entire range.

Centroid shifts

We examined range shifts in 88 species of North American butterfly species by examining changes in abundances across 146 fixed transects between the years 2000 – 2017 (Supplementary Figure 1). Transect monitoring was conducted by citizen scientists that took part in the Ohio Lepidopterists Monitoring Program within a strict range of environmental conditions on a weekly basis (https://www.ohiolepidopterists.org/). We only examined range shifts in species geographic ranges that had over 80 observations across the 18-year test period, and species that had range shift data for at least a 10-year period.

We used abundance-weighted methods to assess each species’ range shifts. We calculated the abundance-weighted geographic centroid of each species each year between 2000 and 2017 within the Ohio Lepidoptera monitoring scheme to create centroid shift vectors for each species. To account for potential sampling error, we bootstrapped species occurrence records by resampling 80% of each species annual records with replacement 10,000 times to estimate each species’ mean geographic centroid each year of the survey. Centroid shift vectors have been shown to be an efficient and accurate way to assess fine scale changes in species distribution, especially over relatively short time scales.

We used linear models to calculate the velocity at which each species shifted their latitudinal and longitudinal centroids across space. The range shift velocities have components of both rate and direction. The slope and standard error associated with each model was extracted to estimate range shift velocity. Degrees latitude were converted to kilometres by multiplying latitude by 84 (the mean degrees to latitude conversion across Ohio's latitudinal span), and degrees longitude were converted to kilometres by multiplying degrees longitude by 111 (same longitudinal conversion across the globe). We calculated velocity of centroid shift as per Huang et al. 2017 equation 1,

sqrt(velocity(lat)squared) + (velocity(long)squared).

Using the slopes of the latitudinal and longitudinal centroid shifts we estimated the start and end points for each species centroid shift within Ohio to visualise range shift vectors. The R package geosphere was used to calculate the distance (km) each species centroid shifted and the direction they shifted (bearing degrees). We plotted the rate and direction of species range shifts using a rose wind diagram from the package clifro.

Climate change velocity 

We calculated the climate change velocities that each butterfly species experienced at each of their annual centroids within Ohio. We did this to examine both how starting climate change velocity impacts range shift inertia, and also to understand what kind of climates butterflies are shifting their ranges towards over time. For the latter question, we calculated the slope of each species’ annual centroid climate change velocities as a function of year. This allowed us to assess the nature of the climate change velocities that species were moving towards (e.g. towards environments that were warming less rapidly and becoming wetter over time) and whether the climates species were moving towards impacted the rate at which they shifted their centroids.  We also extracted the mean climate change velocity species experience throughout their entire geographic distribution using the vocc package. Climate change velocities are calculated as the ratio of temperature change over time (temporal trend), to temperature change over space (spatial trend) based on a 3 x 3 cell neighbourhood, where climate change velocity = spatial trend/temporal trend. We estimated the velocity at which the warmest temperatures in July (hottest month of the year), coldest temperatures in January (coldest month of the year), and mean precipitation received in April (important month for plant growth and butterfly development), changed from 2000 - 2017, matching the period of time that we analyse butterfly range shifts. We also divided Ohio into quadrants and estimated the climate change velocity within each quadrant.

 Niche and dispersal traits

To assess how niche traits might be influencing species range shifts, we extracted occurrence data from the Global Biodiversity Information Facility (GBIF) for each species (list of doi’s for each species range data listed in Supplementary Table 1). We computed species distribution polygons for each species using the rangeBuilder and geosphere  packages. Using the raster package, we estimated the maximum and minimum environmental temperatures species experience within their range polygons at a 30-arc-second resolution to calculate the environmental thermal niche breadth of each species.

As different parts of species ranges have different range shift responses to climate change, we calculated an index of each species geographic range position relative to Ohio. We did this using the minimum and maximum latitudes of each species range polygon to calculate a latitudinal range for each species. We then estimated the quantile of this range that overlapped the latitudinal midpoint of the monitoring scheme. Values range from 0 to 1, with values closer to 1 indicating the monitoring scheme overlaps the species’ range towards its poleward edge and values closer to 0 overlapping toward its equatorward edge.

We estimated dispersal capability based on wingspan (in cm). Minimal and maximal wingspan values from one forewing tip to another were available for the majority of butterfly species, so we tested the predictive ability of both sets of wingspan estimates.

Statistical Analysis 

Statistical analyses were performed in the R version 4.2.1 (R Core Team, 2023). To test the drivers of butterfly range shift responses, we used the ‘rma.mv’ function in the metafor package so that we could account for the error in the rate of centroid shift we estimated when calculating species shift vectors, account for phylogenetic effects, and conduct hypothesis testing.

We used a recently published phylogeny of the butterflies of North America which was constructed using 13 common markers (1 mitochondrial and 12 nuclear) to account for evolutionary history in our models. The phylogeny included 70/88 species that we had centroid shift data for. We examined if there was phylogenetic signal in rate of range shift, wingspan, and thermal niche breadth by calculating Pagel’s 𝜆 using the phytools package. We also examined if there was phylogenetic signal in the direction species shifted their ranges by computing discrete character evolution models using the fitDiscrete function to estimate Pagel’s 𝜆 with the package geiger. We computed equal-rate (ER), symmetric transition (SYM) and all rate different (ARD) evolutionary models with a lambda transformation and then compared them using Akaike Information Criterion (AIC) and likelihood ratio tests.

 We structured our models whereby species rate of range shift was the response variable, and the standard error associated with each species rate of range shift was accounted for. We included local climate change velocity (at the start of each species centroid shift vector) (temperature, precipitation and their interaction) and mean climate change velocity (temperature, precipitation and their interaction) throughout species ranges in the same models to understand how climate change velocity impacts range shift inertia. To avoid introducing multicollinearity into our models, we did not include maximum and minimum temperature velocity in the same model as those models had variance inflation factors (VIFs) over 2.5. Similarly, thermal niche breadth and latitudinal position relative to Ohio were correlated (r2 = 0.64), and models that included both variables had VIFs over 2.5. We elected to include thermal niche breadth in our models as our main hypotheses were focused on this trait; we excluded latitudinal position relative to Ohio to avoid collinearity issues. All models included species and a phylogenetic correlation matrix as random factors to account for shared evolutionary history. We used a maximum likelihood modelling approach in all models. Using model reduction (only variables linked to key hypotheses were included in the full models) and AIC comparison, we examined which models best explained variation in species range shifts. We used likelihood ratio tests for significance testing of the variables in the best-fitting models.

In a separate model, which included the predictor variables from the model with the lowest AIC value, we examined the relationship between rate of range shift and wingspan as a proxy for dispersal capability. Wingspan was analysed separately because wingspan data was only available for a subset of the species in the dataset (n = 66).

We also examined if the rate that species shifted their centroids depended on the direction that butterflies were shifting their centroids, in a model that included rate of centroid shift (and its associated standard error) as a response variable, and direction as a categorical variable (N, S, E, W, NE, NW, SE, SW) (i.e. do butterflies shift their centroids fastest in a certain direction?). Furthermore, to understand if there is a relationship between climate change velocity (local and whole range temperature and precipitation) or wingspan on the direction species shifted their centroids, we conducted separate linear models, where climate change velocity or wingspan were included as continuous response variables and direction was included as a categorical predictor variable (as a categorical variable cannot be used as a response variable).

To understand what kind of climate change velocities butterflies are moving into (i.e. are butterflies moving towards environments that are changing more slowly or more quickly?) and how this impacts the rate at which they shift their centroids, we examined the relationship between rate of centroid shift, and the slope of each species climate change velocity vectors (temperature (maximum and minimum), precipitation and their interaction). We used the same modelling methods and comparisons as described above. Pairing this analysis with the starting centroid vs whole species range climate change velocity analysis allows us to understand how both past and present climate change velocity is shaping species range shift responses.