DOI: 10.5061/dryad.0000000bg
title: Data from: Local climate change velocities and evolutionary history explain multidirectional range shifts in a North American butterfly assemblage

description: 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.
description: 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.
description: # Data from: Local climate change velocities and evolutionary history
 explain multidirectional range shifts in a North American butterfly
 assemblage
 [https://doi.org/10.5061/dryad.0000000bg](https://doi.org/10.5061/dryad.0000000bg) We used citizen science collected data of butterfly occurrence records from the Ohio Lepidopterists Monitoring Program between 2000 - 2017. We estimated species range shifts within Ohio by estimating species annual mean centroids within Ohio and calculating range shift vectors (find detailed methods in the published manuscript).  ## Description of the data and file structure The data file includes columns with the following information:  species_names: scientific names of each species Total.butterfleis.counted: total number of butterflies included in range shift analysis for each species First.year: year of first observation  Last.year: year of last observation year.span: number of years range shift was calculated for sig.shift: whether the range shift was statistically significant (Y/N) shift.velocity: each species total range shift velocity (km/yr) shift.bearing: range shift direction in bearing degree (360).  direction: categorical range shift direction start.centroid.max.temp.vel: maximum temperature velocity at species starting centroid (km/yr) end.centroid.max.temp.vel: maximum temperature velocity at species ending centroid (km/yr) start.centroid.min.temp.vel: minimum temperature velocity at species starting centroid (km/yr) end.centroid.min.temp.vel: minimum temperature velocity at species ending centroid (km/yr) start.centroid.precip.vel: precipitation velocity at species starting centroid (km/yr) end.centroid.precip.vel: precipitation velocity at species ending centroid (km/yr) maxtemp_vel_range_polygon: mean maximum temperature velocity throughout a species range polygon (km/yr) mintemp_vel_range_polygon: mean minimum temperature velocity throughout a species range polygon (km/yr) precip_vel_range_polygon: mean precipitation velocity throughout a species range polygon (km/yr) max.temp.vel.vector.slope: slope of maximum temperature velocity across species annual centroids (km/yr) min.temp.vel.vector.slope: slope of minimum temperature velocity across species annual centroids (km/yr) prec.vel.vec.slope: slope of precipitation velocity across species annual centroids (km/yr) lat.pos: species geographic range position relative to Ohio min.av.temp(bioclim): minimum average temperature of the coldest month species experience throughout their ranges (degrees C) max.av.temp(bioclim): maximum average temperature of the hottest month species experience throughout their ranges (degrees C) Tbreadth: range of thermal environments species experience throughout their ranges (degrees C) ## Sharing/Access information * For the phylogenetic analysis, we used the butterfly phylogeny published by Earl et al. 2021 [Earl, C., Belitz, M. W., Laffan, S. W., Barve, V., Barve, N., Soltis, D. E., Allen, J. M., Soltis, P. S., Mishler, B. D., & Kawahara, A. Y. (2021). Spatial phylogenetics of butterflies in relation to environmental drivers and angiosperm diversity across North America. *Iscience*, *24*(4), 102239.](https://www.zotero.org/google-docs/?MazxmD)