Study region

The study region encompassed 2566 Ordnance Survey 10 km × 10 km grid squares (hectads) covering the British mainland plus any near-shore islands connected to the mainland by the contiguous spread of hectads. During the first recording period (1976–1990), the mean annual temperature was 8.5 °C, increasing to 9.3 °C during the second recording period (2001–2015)36. The level of warming was therefore 0.8 °C (0.03 °C y−1) across the 25-year interval between the midpoints of the two recording periods. Given the latitudinal gradient in mean annual temperature, this equates to a 257–293 km northward shift in isotherms, depending on latitude (Fig. S1). At the median 1976–1990 range margin for our study species, the northward shift in isotherms was 281 km (11.2 km y−1).

Species records

We considered all animal groups represented in the UK Biological Records Centre (BRC, www.brc.ac.uk) and included any group that contained at least five species meeting our inclusion criteria (see sections below). We identified 13 taxonomic groups, all invertebrates, with sufficient data for inclusion (Table S1), including carnivores, herbivores and omnivores, aquatic (freshwater) and terrestrial taxa, groups that disperse in the soil, by walking, by ballooning and by active flight, and spanning orders of magnitude in body mass.

Each group was covered by a formal recording scheme (Table S1). The majority of species records were collected by citizen scientist recorders, before being collated and cleaned by experts in the group/region to filter out possible errors. We retained the taxonomic distinctions and groupings used by these recording schemes (e.g., butterflies and macromoths were treated as separate groups, whereas dragonflies and damselflies were aggregated). Therefore, ‘group effects’ may reflect differences in the recording schemes as well as the effects of taxonomic group per se.

Each biological record represents a unique location × date observation of species presence. We removed records with ambiguous taxonomy (sensu lato, sensu auct, naming multiple species or identified only to genus). Species listed with a sub-species trinomial, with the label sensu stricto, with variants or different morphs were aggregated at the species level. When analysing range shifts, we used all records with at least hectad-level spatial accuracy that could be unambiguously assigned to one of the two recording periods (1976–1990 and 2001–2015). For habitat and climate associations, we used day-specific records accurate to 1-ha resolution across the period 1976–2015 (for the 291 species included in the final analysis, 74% of records had this level of spatial and temporal precision, ranging from 55% in the first recording period to 80% in the second recording period).

Criteria for species inclusion

We selected non-migratory species that reach their northern (cool) range boundaries in southern/lowland Britain. We defined these species as having 90% of their 1976–1990 distribution in the warmest 50% of the study region36 (Fig. S1), with none of these records within 100 km of the north coast. Since these species have historically been concentrated in the warmer half of Britain, it is reasonable to postulate that they might be favoured by climate warming, and that their distributions would be predicted to expand (generally polewards). As non-migrants, the range extensions we document represent the establishment of new populations over multiple generations.

We excluded species classified as non-native, alien-native hybrid, unknown origin, and those that are dependent on non-native species, as defined by the BRC and the GB Non-native Species Information Portal37. We also excluded vagrants and species thought to be extinct from the study region, including species that have been reintroduced following extinction (e.g., Large Blue butterfly Maculinea arion).

This resulted in an eligible species set of 1570 species in 28 groups. Of these, 421 species had sufficient data to calculate range shifts, of which 305 (17 groups) also had sufficient data to calculate habitat associations (criteria detailed below). Excluding groups with fewer than five species resulted in a final dataset of 291 species in 13 groups.

Range-shift calculations

To calculate range shifts, we first controlled for changes in recorder effort over time (1976–1990 to 2001–2015). We restricted distribution data to ‘well-recorded’ hectads, for which at least 10% of the regional species pool for a group was recorded present in both recording periods38 (Fig. S2). For each group × hectad, we defined the regional species pool as the total number of species recorded in the nearest 100 hectads17, using all species in the database for a given taxonomic group (i.e., regardless of the above inclusion criteria).

For all species occupying at least 20 such hectads in both recording periods, we calculated northern (cool) range margins as the mean latitude of the ten-northernmost occupied hectads. We checked that species had sufficient area to expand or retreat from their 1976–1990 range margins. Hence, we excluded any species with fewer than ten recorded (as above) hectads within 100 km to the north of the range margin, and ten such hectads within 100 km to the south of the range margin17. For the remaining species, we defined range shifts as the latitudinal change (km) in range margins between 1976–1990 and 2001–2015. We converted latitudinal changes to annual rates (km y−1) by sampling random dates (10,000 times) from within each of the two 15-year recording periods, dividing latitudinal shifts by each sampled interval, and recording the median rate for each species. When calculating summary rates of range shift across multiple species (e.g., Fig. 1 and Table S2), we did this separately for each sampled time interval and report the 95% confidence interval around the sample median.

Habitat classification

To identify habitat classes, we used a 25 m × 25 m land-cover map for Britain (LCM2007). This map was created by the NERC Centre for Ecology and Hydrology19, using combined summer and winter satellite data (Landsat-TM5, IRS-LISS3, SPOT-4 and SPOT-5 sensors, pixel size of 20–30 m) enhanced with cartographical information (e.g., Ordnance Survey data, soil types, agricultural census boundaries and urban extents). The classification was trained and validated using a large network of habitat surveys and ground reference points, producing an overall accuracy of 83%. Out of 23 habitat classes identified in LCM2007, we discarded one (saltwater), retained 14 as originally mapped, and created four aggregate classes from the remaining eight: ‘heather’ and ‘heather grassland’ became ‘dwarf shrub heath’; ‘supra-littoral rock’ and ‘littoral rock’ became ‘coastal rock’; ‘supra-littoral sediment’ and ‘littoral sediment’ became ‘coastal sediment’; ‘suburban’ and ‘urban’ became ‘built-up and gardens’.

Climate estimates

We estimated climatic conditions corresponding to each 1-ha species record, and the wider range-margin landscape, by spatially downscaling 5 km × 5 km resolution UKCP09 climate grids provided by the UK Met Office36. For each month of the study period (1976–2015), we used universal kriging with linear dependence on elevation to spatially interpolate mean daily minimum (T Min ) and mean daily maximum (T Max ) air temperatures to 1-ha resolution (Fig. S3). Elevation data were from the Ordnance Survey Terrain 50 product, resampled to the 1-ha grid. For each month and climate variable, we constructed spherical and exponential variograms with distance cut-offs of 100 km and retained whichever had the lowest sum of squared errors. Kriging was implemented using the nearest 30 points. Using the same procedure, but with no dependence on elevation, we kriged UKCP09 monthly rainfall and sunshine-hours grids. Sunshine hours were then converted to sunshine fraction (dividing by day length) to estimate the proportion of the day was typically cloud-free in a given month.

Using sunshine fraction and information on topographic position, we adjusted maximum temperatures for solar radiation. First, we used the Solar Area Radiation tool in ESRI ArcMap to calculate solar radiation at 1-ha resolution, once assuming a ‘flat’ surface (SR Flat ) and then again accounting for the influence of slope, aspect and hill shading (SR Topo ). In each case, sky conditions were weighted linearly towards clear skies (transmissivity = 0.7, diffuse fraction = 0.2) or overcast conditions (transmissivity = 0.2, diffuse fraction = 0.7) depending on the sunshine fraction for the corresponding month. We then used the ratio of SR Topo and SR Flat to scale each grid cell’s diurnal temperature range39, as follows:

$${{\rm{T}}}_{{\rm{MaxSR}}}={{\rm{T}}}_{{\rm{Min}}}+{{\rm{SR}}}_{{\rm{Topo}}}/{{\rm{SR}}}_{{\rm{Flat}}}\times ({{\rm{T}}}_{{\rm{Max}}}-{{\rm{T}}}_{{\rm{Min}}}).$$

Mean daily temperatures were calculated as the mean of daily minima and daily maxima. Finally, we derived three biologically-relevant annual climate variables21: minimum winter temperature, degree-days above 5 °C, and the ratio of annual rainfall to potential evapotranspiration40. We defined annual variables over 12-month periods beginning 1 September, because this represents a more biologically-relevant annual cycle in Britain than a calendar year (which begins part-way through winter)35.

Habitat and climate associations

We identified habitat and climate associations using quasibinomial regression of species presence or absence overlaid on the 18 habitat classes (categorical predictor) and the three annual climate variables (continuous predictors)20. The regression equation for each species was used to estimate its probability of occurrence in each habitat class, under the assumption of equal availability of all habitat classes (i.e., as close as is possible to a ‘species characteristic’) and with climate fixed at mean (centred) values. We defined levels of habitat specialisation to be the coefficient of variation across these 18 probabilities18, producing a species’ specialisation score which, for our dataset, ranged from 0.3 (generalist) to 1.9 (specialist). Results are summarised by taxonomic group in Table S3 and reported for individual species in Table S7.

Given the finer grain of the land-cover map (25 m × 25 m) compared with the species data (1 ha), some species records could potentially be associated with multiple habitat classes. In these cases, we assigned the majority habitat class for the 1-ha pixel, and included weights to indicate the proportion of the pixel covered by this habitat41. Climate values were specific to the year of a species observation, except where a species was observed in the same 1-ha pixel in multiple years, in which cases we assigned the mean climate values across those years.

We otherwise took all recorded presences to be ‘true’ for the purposes of modelling and included in the final analysis all species with presence records in at least 200 distinct 1-ha pixels (approximately 20 presence pixels for each parameter estimated by the model; mean records [distinct pixels] per species = 8,167 [2,238], median = 1,390 [694], range = 260–343,040 [201–104,960]). Inferring absence data from presence-only datasets is inherently more difficult. To minimise the number of false absences, we applied the following controls. First, we only included as potential absences those pixels that had been visited by recorders of the same recording scheme (as deduced from records of other species within the same recording scheme), within 50-km of a presence observation for the target species, and excluding landscapes occupied during only one recording period. We did this to account for, respectively, lack of visitation by appropriate recorders, historical barriers to dispersal, and changes in land cover.

Second, we filtered absences according to time of year, given their location, for example to avoid treating late summer data as absences if the target species’ flight period is in spring. We did this by fitting smooth phenology curves to the frequencies of record dates (days of the year) for the target species. To account for spatial variations in phenology, we restricted these records to within 50 km of a particular absence site. Potential absences with record dates in the tails (lower or upper 10% area under the curve) of these location- and species-specific phenology curves were excluded. In the few cases where smoothing splines could not be constructed, we defined the tails of the phenology curve directly from the raw dates (10th and 90th percentiles, correcting for year-breaks where needed).

The remaining absences were from 1-ha pixels visited under the same recording scheme as the target species, in landscapes near where the target species had been recorded in both time periods, and within the appropriate phenological time windows. The absences still varied in reliability, however, because some qualifying pixels had only been visited only once, whereas others had been visited multiple times. Third, therefore, we weighted absence data by the probability of recording the target species if it was present, given the number times (t) the absence pixel was visited:

$$\frac{1}{n}\sum _{s=\mathrm{1..}n}1-{(1-{p}_{s})}^{t}$$

That is, one minus the probability of failing to detect the species on every occasion, where the p s are probabilities of detection across n known presence sites for the target species (calculated as the number of times the species was recorded in site s divided by the number of times s was visited). To account for spatial variations in abundance (and therefore detectability) we calculated p s using data restricted to within 50 km of the absence site.

Conditions at species’ range margins

We obtained spatial estimates of habitat availability by predicting each species’ regression model back on the land-cover map, with climate fixed at mean (centred) values. Habitat availability at the range margin was defined as the mean value across all 25 m × 25 m land-cover pixels in a circular 50-km buffer around the ten-northernmost hectads (or >10 hectads when >1 hectad tied for having the 10th highest latitude) that were used to define the range margin in the first recording period (1976–1990); i.e., landscapes across which the species had potential to expand or retract over time. Habitat availability for individual species ranged from 0.5% (very little of the landscape could be colonised) to 57% of the landscape (ample opportunity for expansion; Tables S4 and S7). Habitat availability exhibited positive skew (Fig. 1; Shapiro-Wilk, W = 0.75, P < 10−20, n = 291 species) and so we log 10 -transformed these values to improve normality (W = 0.997, P = 0.90).

We defined exposure to climate change as the logged ratio of mean range-margin conditions in 1976–1990 versus 2001–2015. That is, we predicted each species’ regression model back on the land-cover map twice, first with annual climate set to its mean condition for the period 1976–1990, and second with climate set to its mean condition for 2001–2015. Since habitat data were constant in the model, the ratio of these predictions was >1 if climatic suitability improved over time, <1 if climatic suitability deteriorated, and =1 if there was no change in climatic suitability. We defined exposure using logged ratios so that there was symmetry about the no-change line; i.e., absolute exposure had the same magnitude when climate was improving as when it was deteriorating. Our usage of the term ‘exposure’ equates to change in climatic suitability for a species and encompasses both sensitivity to climate and the extent to which relevant aspects of climate have changed (cf. IPCC terminology of exposure being independent of sensitivity21).

Models of range shift

We modelled species’ observed range shifts as linear functions of habitat specialisation, log 10 -transformed habitat availability at the range margin, and exposure to climate change at the range margin. We constructed single-predictor models, as well as multivariate models where collinearity among predictors was low (Table S6): habitat specialisation and log 10 -habitat availability were highly correlated (r = −0.53; Fig. 2) and so we did not include both in the same model, whereas exposure to climate change was uncorrelated with both of these predictors (r = −0.02 and r = 0.07, respectively). To account for phylogenetic relatedness and methodological differences in recording between taxonomic groups (i.e., across recording schemes), we used linear mixed-effects models with taxonomic group as the random intercept term (Table S6).

Sensitivity to recording level

We ranked the 13 taxonomic groups by descending geographic coverage of citizen-science recording, defined by the number of hectads where there has been sufficient recording for at least 25% of the regional species richness (considering the nearest 100 hectads) to have been sampled in both recording periods (Table S1). In Fig. 1, we summarise between-species variation separately for: (1) four groups with the high intensity recording; (2) eight groups with lower intensity recording; and (3) macromoths. We plotted macromoths separately because, unlike other groups, moth recording uses attractant methods (light traps at night) so that the areas sampled – and thus habitat associations – are more uncertain.

The proportion of variation in range shift that could be explained was higher for taxonomic groups with higher recording intensity, but the signs of the relationships were consistent (Table S6), demonstrating that the patterns we report are qualitatively robust to recording intensity. In Fig. 4 we systematically varied the threshold of recording coverage, above which species are included in the habitat-climate interaction model. For example, when the recording threshold is low, many groups are eligible for inclusion; when the threshold is high, only the best-recorded groups are included. For consistency in statistical power across different levels of group inclusion, each point in Fig. 4 was generated by averaging over 10,000 randomised draws of 30 species from three qualifying groups.

We conducted statistical and spatial analyses using R version 3.6.1, with packages lme442, lmerTest43, MuMIn44, cAIC422, doParallel45, raster46 and rgeos47.