# Climate-induced phenology shifts linked to range expansions in species with multiple reproductive cycles per year

Callum J Macgregor, Chris D Thomas, David B Roy, Mark A Beaumont, James R Bell, Tom Brereton, Jon R Bridle, Calvin Dytham, Richard Fox, Karl Gotthard, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Callum J Macgregor, Chris D Thomas, David B Roy, Mark A Beaumont, James R Bell, et al. Climate-induced phenology shifts linked to range expansions in species with multiple reproductive cycles per year. Nature Communications, 2019, 10 (1), pp.4455. 10.1038/s41467-019-12479-w. hal-05200392

HAL Id: hal-05200392

https://hal.science/hal-05200392v1

Submitted on 5 Aug 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Climate-induced phenology shifts linked to range expansions in species with multiple reproductive cycles per year

Callum J. Macgregor ID <fnref n="1" />*, Chris D. Thomas ID <fnref n="1" />, David B. Roy ID ², Mark A. Beaumont³, James R. Bell⁴,
Tom Brereton⁵, Jon R. Bridle³, Calvin Dytham ID <fnref n="1" />, Richard Fox ID ⁵, Karl Gotthard ID ⁶, Ary A. Hoffmann⁷,
Geoff Martin⁸, Ian Middlebrook⁵, Sören Nylin⁶, Philip J. Platts ID ⁹, Rita Rasteiro ID ³, Ilik J. Saccheri¹⁰,
Romain Villoutreix¹⁰, Christopher W. Wheat⁶ & Jane K. Hill ID <fnref n="1" />

Advances in phenology (the annual timing of species' life-cycles) in response to climate change are generally viewed as bioindicators of climate change, but have not been considered as predictors of range expansions. Here, we show that phenology advances combine with the number of reproductive cycles per year (voltinism) to shape abundance and distribution trends in 130 species of British Lepidoptera, in response to ~0.5 °C spring-temperature warming between 1995 and 2014. Early adult emergence in warm years resulted in increased within- and between-year population growth for species with multiple reproductive cycles per year (n = 39 multivoltine species). By contrast, early emergence had neutral or negative consequences for species with a single annual reproductive cycle (n = 91 univoltine species), depending on habitat specialisation. We conclude that phenology advances facilitate pole-wards range expansions in species exhibiting plasticity for both phenology and voltinism, but may inhibit expansion by less flexible species.

<footnote n="1">¹Department of Biology, University of York, York YO10 5DD, UK. ²Centre for Ecology and Hydrology, Maclean Building, Benson Lane, Crowmarsh Gifford, Wallingford, Oxfordshire OX10 8BB, UK. ³School of Biological Sciences, University of Bristol, Bristol BS8 1TL, UK. ⁴Rothamsted Insect Survey, Biointeractions and Crop Protection, Rothamsted Research, West Common, Harpenden, Hertfordshire AL5 2JQ, UK. ⁵Butterfly Conservation, Manor Yard, East Lulworth, Wareham, Dorset BH20 5QP, UK. ⁶Department of Zoology, Stockholm University, Stockholm SE-106 91, Sweden. ⁷Department of Zoology and Genetics, Bio21 Institute, University of Melbourne, Parkville 3010 Victoria, Australia. ⁸Department of Life Sciences, Natural History Museum, Cromwell Road, London SW7 5BD, UK. ⁹Department of Environment and Geography, University of York, York YO10 5NG, UK. ¹⁰Institute of Integrative Biology, University of Liverpool, Liverpool L69 7ZB, UK. *email: callumjmacgregor@gmail.com</footnote>
Climate change is resulting in changes in the size, latitudinal range¹,² and elevational extent³ of species' distributions. However, distribution changes are highly variable among species, and rates of polewards expansion often fail to track the climate³⁻⁵. Range expansions are dependent on stable or increasing abundance trends⁶, and hence understanding the effects of climate change on species' abundances is crucial in order to understand variation in range shifts. A potentially important contributing factor is phenological advancement⁷⁻⁹, with many species now undertaking life-cycle events earlier in the year. However, it is unclear whether such phenology advances are beneficial or detrimental for populations of species¹⁰⁻¹⁴.

We used British Lepidoptera to examine this question, because long-term phenological, population and distribution data are all available spanning several decades. We focus on two traits in which Lepidoptera can display phenotypic plasticity (whereby environmental cues directly alter the physical or behavioural phenotype of individuals¹⁵): voltinism and phenology. Lepidoptera include species that are obligately univoltine (i.e., all individuals pass a winter in diapause), and species in which every individual makes a plastic developmental 'decision' whether to undergo diapause or to directly develop, based on environmental cues. Populations of such species may therefore undergo multiple generations per year depending on the length of the local growing season annually, and in cooler regions may be functionally univoltine.

Many Lepidoptera have also advanced their phenology, with adults emerging earlier in recent, warmer years¹⁶ because the growth rate of immature stages increases at warmer temperatures (although photoperiod may regulate phenology in some species¹⁷). Such phenology advances could be either detrimental or beneficial to species, depending on the outcomes of longer or more favourable growing seasons¹⁸⁻²¹ and potential temporal decoupling from host-plants or natural enemies²²⁻²⁴. Overall, it is currently unclear whether phenology advances will result in increases or declines in annual abundance, and whether species with different life histories differ in the consequences of phenology advances.

Here, we show that phenology advances have resulted in increased abundance trends and range expansions in species with multiple reproductive cycles per year (i.e., multivoltine species). Early emergence permits the number of individuals of these species to increase faster in second and subsequent generations within the year, generating positive overall abundance and distribution trends. However, phenology advances do not correlate with abundance trends or range expansions in species with a single annual reproductive cycle (i.e., univoltine species), and are associated with abundance declines in the subset of univoltine species that are also habitat specialists.

## Results

**Interspecific relationships between phenology and demography.** We analysed trends in the phenology, abundance and distribution of 130 species of Lepidoptera (29 butterflies and 101 moths) for which trends could be robustly estimated over a 20-year period (1995–2014) during which mean spring temperatures warmed by ~0.5 °C (Supplementary Fig. 1). We compared functionally univoltine species (n = 91; defined as those that rarely undergo more than one generation per year anywhere in their British range, e.g., Silver-studded Blue *Plebejus argus*; Fig. 1c) with multivoltine species (n = 39; those that regularly undergo two or more generations in part or all of their British range, e.g., Small Blue *Cupido minimus*; Fig. 1d). Both univoltine and multivoltine species significantly advanced their adult emergence dates over the study period, with the first annual emergence peak

for multivoltines (~3 days/decade, range −23.8 to 16.7) advancing significantly faster than univoltines (~1.5 days/decade, range −4.8 to 6.2; Supplementary Table 2).

We found that phenology advances led to positive abundance trends for multivoltine species, but not for univoltine species (likelihood ratio test (LRT), $\chi^2 = 8.23$, d.f. = 1, P = 0.004; Table 1). Multivoltine species showed greater increases in abundance if they had advanced their phenology, but there was no clear relationship between phenology advances and abundance trends among univoltine species (Fig. 2). We found that phenology advances did not directly correlate with change in the distribution size (LRT, $\chi^2 = 0.07$, d.f. = 1, P = 0.792) or change in range margin latitude (LRT, $\chi^2 = 1.29$, d.f. = 1, P = 0.256). However, abundance trends were themselves significantly, positively related to trends in both distribution size (LRT, $\chi^2 = 52.3$, d.f. = 1, P < 0.001) and range margin (LRT, $\chi^2 = 8.82$, d.f. = 1, P = 0.003) for all species, regardless of voltinism (Supplementary Fig. 2). To test the indirect relationship between distribution and phenology moderated by abundance, we predicted species' abundance trends from our models of the relationship between phenology advances, voltinism and abundance trends (Table 1), yielding an estimate of the specific component of abundance change that was driven by phenology advances. We found that these model-predicted abundance trends were significantly related to trends in range margin (LRT, $\chi^2 = 5.16$, d.f. = 1, P = 0.023), but marginally not distribution size (LRT, $\chi^2 = 3.17$, d.f. = 1, P = 0.075). Hence, we conclude that climate-linked shifts in range margin latitudes are indirectly driven by phenology advances, mediated by effects of abundance (Fig. 1). These results were robust to phylogeny, to the designation of some species (n = 36) whose voltinism patterns were hard to categorise because they have both univoltine and bivoltine populations in Britain and/or mainland Europe, and to our selection of a relatively short 20-year study period (Tables S3, 4).

**Intraspecific relationships between phenology and demography.** To understand whether links between phenology advances and abundance trends were causally related both between and within species, or solely correlated at species-level, we calculated trends in phenology and abundance independently for every population (i.e., recording site) in our data set (n = 3–104 populations per species) and assessed the intraspecific relationships between the two variables at the population level. These analyses confirmed our previous findings, revealing that multivoltine species showed greater increases in abundance in populations that had advanced their phenology, but no such effect in populations of univoltine species (Table 2). Among multivoltine species, 8/39 (20.5%) species showed individually significant, positive population-level relationships between phenology advances and abundance trends (eight times higher than the two-tailed chance expectation), no species displayed a significant negative relationship, and the average relationship across all 39 multivoltine species was significantly positive (LRT, $\chi^2 = 57.50$, d.f. = 1, P < 0.001; Fig. 2d). These positive population-level relationships were evident in both the subset of multivoltine species that are, on average, increasing nationally, as well as the multivoltine species that are declining nationally (Supplementary Fig. 3), suggesting that phenology shifts could be locally adaptive by limiting population-level rates of decline in species that are declining nationally. By contrast, among univoltine species, only 3/91 (3.3%, c.f. null expectation of 2.5%) species displayed a significant positive relationship between phenology advance and abundance trends at the population level, with 3/91 (3.3%) showing a negative relationship, and the average relationship across all 91
<figure><img src="image_4.png" /></figure>

Low life-cycle plasticity
(univoltine species)

<figure><img src="image_5.png" /><figcaption>a</figcaption></figure>

High life-cycle plasticity
(multivoltine species)

<figure><img src="image_6.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Cupidoid minimus</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>f</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>g</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>h</figcaption></figure>

<figure><img src="image_13.png" /></figure>

univoltine species was not significantly different from zero (Fig. 2c).

**Annual phenological variation.** To gain insight into why multivoltine species benefitted from phenology advances, but

univoltine species did not, we examined the effects of annual variation in emergence dates on a sequence of Lepidopteran life-cycle events (Table 1). Both univoltine and multivoltine species emerged significantly earlier in years when spring temperatures were warmer (Fig. 3a), but with a significantly larger effect in
**Fig. 1** Effects of phenology advances on abundance and distribution trends depend on voltinism, illustrated with example species. The flowchart describes the main conclusions of this study: climate-driven phenology advances have a positive effect on abundance trends in multivoltine species, but a neutral or negative effect on abundance trends in univoltine species (depending on habitat specialisation). In turn, abundance trends have a positive effect on trends in distribution size and northern range margin, regardless of voltinism. Trends in emergence date, abundance, distribution size and northern range margin are depicted for two butterflies: Silver-studded Blue *Plebejus argus*, and Small Blue *Cupido minimus*. **a, b** *P. argus* has retracted in distribution size (−0.2%/yr) and range margin (−1.8 km/yr), whereas *C. minimus* has expanded in distribution size (1.3%/yr) and range margin (7.4 km/yr). Distribution size depicted as number occupied hectads in 1995–2014 (black circles), 1995–2004 only (orange circles) and 2005–2014 only (blue circles). Range margins are depicted for 1995 (orange lines) and 2014 (dashed black lines). **c, d** Voltinism of univoltine *P. argus* and multivoltine *C. minimus* shown by observed abundance on transect counts across all sites and years (grey circles; counts > 100 omitted) and GAM-fitted curves (blue lines). **e–h** Both *P. argus* and *C. minimus* have advanced their phenology (0.25 days/yr and 0.46 days/yr, respectively); *P. argus* has declined in abundance (−5.5%/yr), but *C. minimus* has increased in abundance (3.8%/yr). Observed peak day of first-generation emergence (**e, f**) and mean abundance per recording event (**g, h**) in each year is shown (grey circles), and points at the same site are connected (grey lines); overall trend across the duration of the study period (Supplementary Data 1) is shown (black dashed line) with 95% confidence intervals (grey shading)

multivoltine than univoltine species. Earlier emergence in multivoltine species was associated with greater population growth between the first and second generations (Fig. 3b), and consequently, earlier emergence by multivoltine species led to increased abundance, in both years *t* and *t* + 1 (Fig. 3c, d). However, earlier emergence in univoltine species was significantly associated with reduced abundance in year *t* + 1 (Fig. 3).

**Role of habitat specialisation.** Habitat availability can be an important predictor of range expansion<sup>6,25</sup>, and we found that including habitat specialisation in our statistical models revealed further distinctions between habitat specialist species (*n* = 21) and wider-countryside generalist species (*n* = 109). Most notably, phenology advances led to abundance declines among univoltine habitat specialists, but there was no relationship among univoltine wider-countryside generalists (Table 1). After refining our models by including habitat specialisation, species’ model-predicted abundance trends were significantly related to trends in distribution size (LRT, $\chi^2$ = 14.49, d.f. = 1, *P* < 0.001) as well as range-margin latitude (LRT, $\chi^2$ = 5.90, d.f. = 1, *P* = 0.015).

Discussion

Our results demonstrate that positive demographic responses to climate change are only evident in the subset of species which are multivoltine and have advanced their phenology, thereby showing plasticity in both phenology and voltinism. This combination provides a pathway by which benefits can be gained from earlier emergence in warmer springs, yielding increasing abundance in the second annual generation associated with these phenology advances. Such benefits are not experienced by the subset of multivoltine species that have not advanced their phenology, whilst univoltine species are constrained to develop through only one generation per year (by innate factors or by climate, in the case of species that are functionally univoltine under cool British conditions, even if they have the potential to be multivoltine elsewhere).

Phenology advances are associated with declines in abundance among univoltine habitat specialists, which might experience direct costs associated either with phenology advances themselves, or with warmer winter or spring temperatures (e.g., extended periods of delayed sexual maturity<sup>26</sup>, diapause or larval<sup>27</sup> or adult aestivation<sup>28</sup>), before suitable conditions for emergence/reproduction arise the following year. These factors may potentially lead to greater reproductive success in years of later emergence (i.e., phenology delays; Fig. 3c, d), even though abundance was higher in individual years of earlier emergence (Table 1). In particular, univoltine habitat specialists (whose host–plant niche is often narrow) may experience phenological mismatches with host–plants<sup>29</sup>, from which generalist species

may be buffered. Host–plants may also be advancing their phenology<sup>8</sup>, so it is possible that such univoltine species might have declined even more without phenology advances. As with plasticity in both phenology and voltinism, spatiotemporal variation in habitat and host–plant associations (and hence specialisation) may also include an element of plasticity, whereby individuals can make behavioural decisions to occupy favourable habitats under specific environmental conditions<sup>30</sup>; but genotypic diversity among populations and individuals may also contribute to habitat and host–plant selection.

Despite this, we found no significant differences between univoltine and multivoltine species in their overall abundance or distribution trends, potentially because emergence dates did not advance for all multivoltine species, or for all populations of species that have advanced their phenology overall (Supplementary Fig. 3). Identifying the factors that drive or constrain phenology advances will therefore be important. These factors might include local adaptation to photoperiod signals<sup>31,32</sup>, physiological barriers limiting increases in development rate<sup>33</sup>, or availability of suitable habitat or microclimatic conditions to otherwise mitigate the effects of climate change through behavioural responses<sup>30</sup>. In particular, understanding drivers of detrimental phenology advances in populations of univoltine habitat specialists may be important for their conservation, particularly as this group includes many UK conservation-priority species (e.g., High Brown Fritillary *Argynnis adippe*). Our findings show that the effects of climate-driven phenology advances have the potential to be detrimental in univoltine habitat specialists, but some of these species might benefit from gaining a second generation in Britain under future climate change (those which are functionally univoltine in Britain but have the capacity to be multivoltine; e.g., *P. argus*, Fig. 1). This may be more likely to occur if emergence dates continue to advance in much warmer years. Given that there is not a consistent outcome of phenology advances among species, strategies for conservation management under climate change should employ approaches that generate local conditions for a diverse range of phenological strategies across species at different trophic levels, such as approaches that maximise habitat, microclimate and host–plant heterogeneity<sup>34,35</sup>.

In conclusion, our study shows that range expansions in response to climate change<sup>3,5</sup> are influenced by phenology advances, through their effect on population abundance. However, the nature of the relationship between phenology advances and abundance depends on life-cycle plasticity in voltinism. Species with multiple reproductive cycles per year may be able to capitalise on warmer springs by advancing their phenology, thus increasing the total number of reproductive cycles per year<sup>19,20</sup> and/or increasing reproductive success within each cycle<sup>18</sup>, with consequent population growth and expanding distributions. By contrast, univoltine habitat specialists experience apparent costs
<figure><img src="image_14.png" /></figure>

**Fig. 2** Phenology advances and voltinism drive abundance trends in multivoltine species. Advancing phenology correlates with increasing abundance at species- and population-level for multivoltine species (n = 39), but not univoltine species (n = 91). Lines depict model-predicted relationships between phenology and abundance trends, from generalised linear mixed-effects models. **a, b** Points show changes in phenology and abundance over the study period (1995–2014) at species-level for univoltine and multivoltine species, respectively. Point colour indicates taxonomic group (butterflies: filled, moths: open). Lines show significant (*P* < 0.05) relationships. **c, d** Points show changes in phenology and abundance over the study period at the population level. Grey lines show relationships calculated independently for each species; lines are solid if the relationship is significantly different to zero (*P* < 0.05), or otherwise are dotted. Solid black lines indicate the overall relationship across species, and are plotted only if significant (*P* < 0.05). Among univoltine species (**c**), 3/91 (3.3%) show a significant positive relationship, and 3/91 (3.3%) a significant negative relationship, between phenology and abundance change; the average relationship across species is not significant (Table 2). Among multivoltine species (**d**), 8/39 (20.5%) show a significant positive relationship between phenology and abundance change, and none show a significant negative relationship; the average relationship across species is also significantly positive (Table 2)

when they advance their phenology, resulting in population declines and retracting distributions, although these costs are not evident in univoltine generalists. These variable demographic consequences of phenology advances may help to explain why some species' distributions have not expanded quickly enough to track temperature changes³,⁵.

## Methods

**Data sets.** We used data obtained by four recording schemes to assess changes in phenology, abundance, distribution size and latitude of the northern range margin over a 20-year period (1995–2014). Specifically, we used data from two population-monitoring schemes that contained abundance records, with high spatial and temporal resolution over many years for fixed sampling locations, to measure phenology and abundance (butterflies: the United Kingdom Butterfly Monitoring Scheme (UKBMS); moths: the Rothamsted Insect Survey (RIS) Light-Trap Network). We also used data from two distribution recording schemes that contained annual presence records summarised at hectad (10 × 10 km) level for the whole of Great Britain to measure distribution size and range margin (butterflies: Butterflies for the New Millennium (BNM); moths: the National Moth Recording Scheme (NMRS)).

In the UKBMS, data are collected annually over a 26-week period (1st April–29th September). Weekly transects are walked along a fixed route, following a standard method³⁶ to count the abundance of each species present. In the RIS,

night-flying and crepuscular moths are attracted to a 200 W tungsten bulb installed within a standard light-trap design, operated in the same location on every night of the year between dusk and dawn³⁷. Sampled moths are collected daily or every few days, and the abundance of each species counted. Thereby, both recording schemes generate abundance data for a fixed site, with a temporal resolution of 1 week or better, over a long period of time (in many cases continuously for 2–4 decades). This allows for reliable estimation of changes in site-level abundance and phenology over time.

In both the BNM and the NMRS, data are contributed with high spatial resolution by volunteer recorders as a form of citizen science, and summarised to produce annual distribution maps at hectad resolution. The BNM was established in 1995³⁸ but builds upon a previous atlas project³⁹, whilst the NMRS officially commenced in 2007; both recording schemes include historical records dating back to the 17th and 18th centuries, respectively. Both schemes comprise mainly records of adult Lepidoptera, either observed during the daytime or captured in light traps, but also include other recording methods (e.g., pheromone lures) and records of immature life-stages. The annual number of the hectad-level species presences recorded by the BNM has remained roughly stable in each year since its commencement in 1995 (with fewer records from earlier years), despite growth in the total number of records submitted to the scheme. Both the number of hectad-level species presences recorded by the NMRS, and the number of records submitted, continue to grow⁴⁰. Not all data sets had been updated beyond 2014 at the time of analysis; therefore, we selected the 20-year study period 1995–2014 in order to contain the maximum informational content across the four data sets.
<table><caption>Table 2 Statistical tests of intraspecific relationships between phenology and abundance change</caption><thead><tr><th>Dependent variable</th><th>Independent variable</th><th>Data subset</th><th>Interacting covariates</th><th>n</th><th>AIC</th><th>Marginal R<sup>2</sup></th><th>Effect size (s.e.)</th><th>X<sup>2</sup> (P)</th></tr></thead><tbody><tr><td rowspan="10">Change in abundance</td><td rowspan="10">Change in emergence date</td><td rowspan="2">Full data set</td><td>Voltinism</td><td>1677</td><td>-5255.4</td><td>0.028</td><td>-</td><td><strong>9.18</strong> (<strong>0.002</strong>)</td></tr><tr><td>Voltinism*</td><td>1677</td><td>-5261.4</td><td>0.032</td><td>-</td><td>0.17 (0.679)</td></tr><tr><td rowspan="2">Univoltine species</td><td>Class</td><td>-</td><td>-3182.7</td><td>0.000</td><td>-0.00 (0.00)</td><td>0.02 (0.886)</td></tr><tr><td>Class</td><td>-</td><td>-3179.4</td><td>0.002</td><td>-</td><td>0.81 (0.368)</td></tr><tr><td>Univoltine habitat specialists</td><td>-</td><td>99</td><td>-254.9</td><td>0.007</td><td>0.01 (0.02)</td><td>0.59 (0.441)</td></tr><tr><td>Univoltine wider-countryside generalists</td><td>-</td><td>939</td><td>-2933.2</td><td>0.0001</td><td>-0.00 (0.00)</td><td>0.13 (0.723)</td></tr><tr><td rowspan="2">Multivoltine species</td><td>Class</td><td>639</td><td>-2029.7</td><td>0.074</td><td>0.01 (0.00)</td><td><strong>57.50</strong> (<strong>&lt;0.001</strong>)</td></tr><tr><td>Class</td><td>639</td><td>-2086.9</td><td>0.084</td><td>-</td><td>0.006 (0.938)</td></tr><tr><td>Univoltine habitat specialists</td><td>-</td><td>14</td><td>-32.46</td><td>0.002</td><td>0.01 (0.04)</td><td>0.06 (0.800)</td></tr><tr><td>Univoltine wider-countryside generalists</td><td>-</td><td>625</td><td>-1995.7</td><td>0.078</td><td>0.01 (0.00)</td><td><strong>58.63</strong> (<strong>&lt;0.001</strong>)</td></tr><tr><td>Multivoltine species (increasing abundance)</td><td>-</td><td>295</td><td>-903.9</td><td>0.068</td><td>0.01 (0.00)</td><td><strong>21.46</strong> (<strong>&lt;0.001</strong>)</td></tr><tr><td>Multivoltine species (declining abundance)</td><td>-</td><td>344</td><td>-1162.7</td><td>0.091</td><td>0.01 (0.00)</td><td><strong>33.80</strong> (<strong>&lt;0.001</strong>)</td></tr></tbody></table>

An overall model was constructed in each case with species as a random effect, and with voltinism as an interacting covariate if indicated. Significance of each model was tested using a likelihood ratio test, and tests that had statistical significance (P < 0.05) are indicated in bold

**Data selection.** To obtain consistent estimates of variables across data sets, unbiased by increased recording in later years, we restricted each data set according to uniform criteria for both butterflies and moths. We grouped subspecies at the specific level by reference to a recent checklist of British Lepidoptera<sup>41</sup>, and treated species complexes as a single taxonomic entity equivalent to one species (the only such aggregate included our final data set was Common/Lesser Common Rustic *Mesapamea secalis/didyma*). We initially excluded from the study: (i) species that are obligatory migrants, or for which a substantial proportion of records represent immigrant individuals; and (ii) species for which new methods of recording have been developed within the study period (e.g., Sesiidae, now mainly recorded using pheromone lures).

For the population-monitoring schemes, we first restricted each data set to include a population (defined as one species at one transect/trap location) in each year only if (i) there were at least ten recording events in that year during which any species was recorded (even if the focal species was not) and (ii) the focal species was itself recorded during at least three of those recording events. For all remaining combinations of population × year, we fitted a generalised additive model (GAM) to all abundance records (including zeroes), with a Poisson error distribution and using a restricted maximum likelihood approach to estimation of smoothing. We followed a series of logical steps (Supplementary Fig. 7) to exclude GAMs which were deemed not to have fitted successfully; GAMs were discarded if their predicted abundance on 1st January was >1, or failed to reach at least one peak (defined as a day on which model-predicted abundance was greater than both the preceding and following days) before 31st December. These rules excluded populations in each year from which first-generation individuals were recorded on the first or last day of recording (UKBMS: April–September; RIS: January–December), preventing reliable estimation of phenology. We then further restricted the data set to include only populations that had successfully fitted GAMs for at least (i) 15 years of the 20-year study period, and (ii) 1 year in the period 1980–1990 (to exclude sites that were recently colonised at the start of the study period, potentially influencing abundance trends<sup>42</sup>). Finally, from the remaining populations, we included only species for which (i) at least three populations matched the criteria above, and (ii) records existed in each of the 20 years of the study period from at least one population, even if no single population had been recorded for 20 years. For some butterflies which may be active before the commencement of UKBMS monitoring (e.g., Peacock *Aglais io*), a sufficient proportion of GAMs fitted unsuccessfully that too little data remained for the species' inclusion in the study, despite being common and widespread; others (e.g., Orange-tip *Anthocharis cardamines*) are represented by only a few populations. It is possible that these early-emerging species may have experienced some of the largest phenology advances<sup>22</sup>.

For the distribution recording schemes, we first restricted each data set separately to include only hectads that were heavily recorded, following previous studies<sup>1,2</sup>, in order to be confident that species not recorded in a hectad were truly absent. Specifically, we first excluded hectads unless they had a presence record (of any species) in both the first and second halves of the study period (i.e., 1995–2004 and 2005–2014, respectively). For each remaining hectad, we calculated annual species richness of the hectad itself and of the 100 nearest neighbouring hectads combined (i.e., the surrounding region), and from these, the annual percentage of

regional species richness that had been recorded in each hectad. We excluded all hectads for which the median annual percentage of regional species richness recorded (across all 20 years) was <25%. This left 1639 heavily recorded hectads in the BNM data set and 475 heavily recorded hectads in the NMRS data set (Supplementary Table 1). Using distributions from within the remaining hectads, we excluded: (i) species which had been recorded in <20 heavily recorded hectads across the full study period (e.g., Lulworth Skipper *Thymelicus acteon*, Dark Bordered Beauty *Epione vespertaria*), as these distributions were too small to reliably estimate change in distribution size; and (ii) species for which the mean elevation of all recorded hectads was >200 m (using elevation data from Farr et al.<sup>43</sup>; e.g., Scotch Argus *Erebia aethiops*, Scotch Annulet *Gnophos obfuscata*), as responses to climate change in upland species might involve elevational shifts rather than changes in distribution or abundance<sup>3</sup>. Finally, we assessed which species reached their range margin >100 km south of the northernmost point of mainland Great Britain (latitude 58° 38' 14" N), following Hickling et al.<sup>1</sup>. We excluded these northerly or ubiquitous species (e.g., Meadow Brown *Maniola jurtina*, Dark Arches *Apamea monoglypha*) from a subset of data for the specific analysis of range margin trends, but retained them in the main data set.

The remaining, final data set contained 130 species for which we had retained reliable data from both population and distribution recording schemes, including 29 butterflies and 101 moths, which represents ~50 and 15%, respectively of all resident British butterfly and macro-moth species (moths are probably more likely to go unrecorded at a site in any given year, despite continuous presence, leading to a lower proportion of populations meeting the requirement for having been recorded in 15/20 years). Of these species, 12 butterflies (41%) were habitat specialists, but only 9 moths (9%), probably because UKBMS transects are more likely to be established on priority habitats occupied by habitat specialists (e.g., calcareous grassland) than RIS traps. For the population-monitoring schemes, our data set comprised 425,087 abundance records of 3,484,983 individual Lepidoptera, spanning 1685 populations at 141 different sites (Supplementary Table 1). From the distribution-recording schemes, it comprised 913,037 hectad-level presence records from heavily recorded hectads.

**Generation of variables.** We generated two categorical variables to describe each species' life-cycle plasticity, by reference to commonly used identification resources<sup>44–48</sup>. First, we described species' functional voltinism: species were classified as univoltine if they rarely or never undergo more than one generation per year anywhere in any part of their range within Great Britain, even if they have the capacity to do so elsewhere in their global range (e.g., Silver-studded Blue *Plebejus argus*, Fig. 1c), or otherwise multivoltine if they regularly undergo a substantial second generation in any part of their British range (e.g., Small Blue *Cupido minimus*, Fig. 1d). In total, we categorised 91 species as univoltine and 39 species as multivoltine. We additionally recorded whether this categorisation was considered to be representative of all populations in all years. Second, we described species' habitat specialisation: species were classified either as habitat specialist or wider-countryside generalist. Butterfly assignments were drawn directly from an established set of habitat specialisation classifications<sup>48</sup>, and moths were assigned by expert opinion, using the same criteria (Supplementary Table 5). In total, we
<figure><img src="image_15.png" /></figure>

**Fig. 3** Species' responses to annual variation in spring temperature depend on voltinism. Lines depict model-predicted relationships ± 95% confidence interval, from generalised linear mixed-effects models fitted to annual population-level estimates of emergence date and abundance for all univoltine and multivoltine species. The relative density of underlying data points is represented by contour lines. Colour indicates voltinism (blue: univoltine species, orange: multivoltine species). **a** Average peak day of first-generation emergence is earlier in years with warmer springs (measured as the accumulated growing degree days above 5 °C (GDD5) between 1st March and 31st May) for both univoltine and multivoltine species, and the effect is significantly stronger for multivoltine species. **b** In multivoltine species, abundance in second and subsequent generations is proportionally larger compared with the first generation (as indicated by a larger intergenerational abundance ratio) in years when peak day of first-generation emergence was earlier. **c** Abundance of multivoltine species (measured as the mean number of individuals recorded per transect/trap) is greater in years when peak day of first-generation emergence was earlier, but there is no relationship for univoltine species. **d** Abundance of multivoltine species is greater when peak day of first generation in the previous year was earlier, but abundance of univoltine species is greater when peak day of first-generation emergence in the previous year was later

categorised 109 species as wider-countryside generalists and 21 species as habitat specialists. The majority of habitat-specialist species were butterflies, despite most species in the overall data set being moths (Supplementary Table 6), reflecting a greater tendency for UKBMS transects to be established in protected areas (with associated habitat specialists)⁴⁹ than for RIS traps.

For all species, we generated four annual variables: first-generation emergence date (phenology), abundance, occupied distribution size and range margin latitude. Phenology and abundance were calculated separately for every population, using population-monitoring scheme data. Single national values for distribution and range margin were calculated in each year, using distribution recording scheme data.

Abundance of each population was the total number of individuals recorded in each year, divided by the number of recording events (transects walked or trap samples collected) in that year, and therefore represented mean abundance per recording event. To estimate the annual phenology of each population, we used the GAM fitted to abundance data (as described above). We used a series of logical steps (Supplementary Fig. 4) to identify the most plausible date for the first peak in abundance; phenology therefore refers to the emergence of the first generation in each year, regardless of each species' voltinism. We used this approach to estimating phenology because it is more robust to the influence of variation in abundance than other approaches (e.g., first appearance date)⁵⁰.

For multivoltine species, we additionally estimated the ratio between abundance in the first generation and all subsequent generations for multivoltine species (intergenerational abundance ratio). We used logical steps again (Supplementary Fig. 4) to identify the most plausible date for the trough of minimum abundance between the first and second generations, and calculated the ratio between the sum of daily abundances (predicted from the GAM) before and after this trough.

To calculate each species' annual distribution size, we calculated the number of heavily recorded hectads in which the focal species was recorded, and the total annual number of heavily recorded hectads in which any species was recorded. From these, we calculated the percentage of the maximum possible distribution size that was occupied by the focal species in each year (distribution). This accounts for an increase in the number of heavily recorded hectads that were recorded in later years of the study period. To calculate the annual latitude of each species' range margin, we identified the ten most northerly occupied hectads (including all hectads that were tied for 10th place) in each year, and calculated the mean northing of these hectads.

We then calculated rates of change over 20 years in phenology, abundance, distribution and range margin for each species (Supplementary Data 1). Change in distribution was calculated as the slope of a linear regression between distribution and year, and was therefore the annual change in the percentage of hectads that were occupied. Likewise, change in range margin was calculated as the slope of a linear regression between range margin and year, and was the annual northwards advance in the latitude of the range margin, in km (a negative value indicated a southwards retraction). Change in abundance was calculated as the slope of a generalised linear mixed-effects model (GLMM) between the logarithm of mean abundance (per recording event) and year, with site as a random effect and a Gaussian error distribution, and was therefore the annual change in abundance as the logarithm of the odds ratio. Finally, change in phenology was calculated as the slope of a GLMM between phenology and year (with the same structure as above), and was therefore the annual change in phenology in days. We reversed the sign of this slope, so that a positive number indicated an advance in phenology (emerging earlier in the year). For the variables generated from population-monitoring scheme data (phenology advance and abundance trend), we additionally calculated
the rate of change separately for each population. These were calculated as above, except that they were the slope of a linear regression rather than a GLMM. For change in abundance and change in phenology, we also calculated rates of change (as above) over the full time period of available data for each species (31–44 years per species, between 1973 and 2017).

We tested the relevance of our species-level trends, based on a subset of data-rich populations, to national trends, using abundance as a case study. Our estimated abundance trends were significantly correlated to long-term national abundance trends from the UKBMS⁵¹ (1976–2016; *F*-test, adjusted *R*² = 30.0%, *F* = 12.97, *P* = 0.001) and RIS⁵² (1968–2016; *F*-test, adjusted *R*² = 31.2%, *F* = 46.29, *P* < 0.001). For the UKBMS, we also calculated national trends for the study period only; these correlated even more strongly with our estimated trends (1995–2014; *F*-test, adjusted *R*² = 49.3%, *F* = 28.25, *P* < 0.001).

Finally, for each of the 141 sites from which we had population-monitoring scheme data included in the final data set, we calculated annual spring temperatures, as the number of growing degree days above a 5°C threshold (GDD5) from 1st March to 31st May inclusive, for the 5 × 5 km grid square containing the site centroid, using gridded data from the UK Meteorological Office⁵³.

**Statistical analysis.** We used GLMMs to test relationships between change in phenology, abundance, distribution and range margin. In each case, we initially constructed two models on the full data set (n = 130 species), testing the fixed effects, respectively, of a two-way interaction between the independent variable and voltinism, and of a three-way interaction between the independent variable, voltinism and habitat specialisation class. We used Gaussian error distributions because our dependent variables were all approximately normally distributed (Supplementary Figs. 5–7), and included taxon group (butterfly or moth) as a random effect (allowing random intercepts). We tested significance of fixed effects in each model using likelihood ratio tests; where interaction terms were non-significant, we retested models with them removed and their constituent parts included, first as two-way interaction terms and if still non-significant, as single main effects. If the final model contained a significant interaction term, we split the data set into subcategories as indicated by the interaction term, and tested whether the relationship between independent and dependent variables was significantly different to zero separately for each subcategory except multivoltine habitat specialist species, because the subset of data for this category was too small (n = 3 species).

Using this approach, we first tested the interspecific effect of change in phenology on all three main dependent variables, calculated at the species level: change in abundance, distribution and range margin. We additionally repeated these analyses using two subsets of data, first excluding species for which the voltinism classification might not be representative of some populations, and second only including such species. Next, we repeated the initial analyses, using a three-level categorical variable to describe voltinism (obligate univoltine, functionally univoltine and multivoltine), where species which have the capacity to be multivoltine but are functionally univoltine throughout Britain were assigned to a separate category from species which are univoltine throughout their global range. Finally, we retested the relationship between change in phenology and change in abundance, using trends in each variable calculated over the full time period of available data.

Next, we tested the direct effects of change in abundance on change in distribution and range margin, because earlier studies suggest that distribution expansions are dependent on stable or positive abundance trends⁶. We also hypothesised that the effects of change in emergence date might indirectly explain change in distribution, mediated by change in abundance, so we used our earlier models to predict the expected change in abundance of each species, based on (i) its voltinism and observed change in phenology, and (ii) voltinism, habitat specialisation, and observed change in phenology. We tested the relationship between these model-predicted changes in abundance, and change in distribution and range margin.

To check for the possible influence of phylogenetic relatedness on our results, we retested these interspecific relationships using phylogenetic generalised least squares (PGLS) models. For this purpose, we constructed a phylogeny of all 130 study species (Supplementary Data 2), using the marker *cytochrome c oxidase subunit I* (COI). We visually confirmed that (i) relationships between families within our phylogeny broadly matched a recent published phylogeny of the Lepidoptera⁵⁴ and (ii) congeneric species were always grouped in monophyletic taxa within our phylogeny.

Thirdly, we tested the intraspecific effect of change in phenology on change in abundance, calculated at the population level, using species as a random effect in place of the taxon group. Finally, we used the annual estimates of spring temperature, phenology, mean abundance per recording event and intergenerational abundance ratio to conduct several further tests, using the same approach as above except with species and year as crossed random effects. When analysing mean abundance and the intergenerational abundance ratio as dependent variables, we used the logarithm of each variable. Specifically, we analysed: (i) the effect of spring temperature upon phenology; and the effect of phenology upon (ii) mean abundance in the same year and (iii) in the following year; and (iv) the intergenerational abundance ratio (for multivoltine species only, because this variable could not be estimated for univoltine species).

All statistical analyses were conducted in R version 3.5.0⁵⁵, except construction of the phylogeny for PGLS, which was conducted in Geneious version 11.1.4⁵⁶. We used the following R packages: mgcv⁵⁷ to fit GAMs, lme4⁵⁸ to construct and test GLMMs, caper⁵⁹ to construct and test PGLSs, ggplot2⁶⁰ to prepare figures and slightly⁶¹ to plot maps of the UK shown in Fig. 1.

**Reporting summary.** Further information on research design is available in the Nature Research Reporting Summary linked to this article.

## Data availability

Data sets were obtained, respectively, from the UKBMS, Rothamsted Research (RIS) and Butterfly Conservation (BNM and NMRS), and may be requested from the same sources.

## Code availability

All R scripts, from initial processing of data sets to final analyses, are archived online at Zenodo.

Received: 8 April 2019; Accepted: 9 September 2019;

Published online: 24 October 2019

## References

1. Hickling, R., Roy, D. B., Hill, J. K., Fox, R. & Thomas, C. D. The distributions of a wide range of taxonomic groups are expanding polewards. *Glob. Chang. Biol.* **12**, 450–455 (2006).

2. Mason, S. C. et al. Geographical range margins of many taxonomic groups continue to shift polewards. *Biol. J. Linn. Soc. Lond.* **115**, 586–597 (2015).

3. Chen, I.-C., Hill, J. K., Ohlemüller, R., Roy, D. B. & Thomas, C. D. Rapid range shifts of species associated with high levels of climate warming. *Science* **333**, 1024–1026 (2011).

4. Mair, L. et al. Temporal variation in responses of species to four decades of climate warming. *Glob. Chang. Biol.* **18**, 2439–2447 (2012).

5. Devictor, V. et al. Differences in the climatic debts of birds and butterflies at a continental scale. *Nat. Clim. Chang.* **2**, 121 (2012).

6. Mair, L. et al. Abundance changes and habitat availability drive species' responses to climate change. *Nat. Clim. Chang.* **4**, 127–131 (2014).

7. Cohen, J. M., Lajeunesse, M. J. & Rohr, J. R. A global synthesis of animal phenological responses to climate change. *Nat. Clim. Chang.* **8**, 224–228 (2018).

8. Walther, G.-R. et al. Ecological responses to recent climate change. *Nature* **416**, 389–395 (2002).

9. Kauserud, H. et al. Climate change and spring-fruiting fungi. *Proc. Biol. Sci.* **277**, 1169–1177 (2010).

10. Thackeray, S. J. et al. Phenological sensitivity to climate across taxa and trophic levels. *Nature* **535**, 241–245 (2016).

11. Both, C., van Asch, M., Bijlsma, R. G., van den Burg, A. B. & Visser, M. E. Climate change and unequal phenological changes across four trophic levels: constraints or adaptations? *J. Anim. Ecol.* **78**, 73–83 (2009).

12. Franks, S. E. et al. The sensitivity of breeding songbirds to changes in seasonal timing is linked to population change but cannot be directly attributed to the effects of trophic asynchrony on productivity. *Glob. Chang. Biol.* **24**, 957–971 (2018).

13. Reed, T. E., Grotan, V., Jenouvrier, S., Sæther, B.-E. & Visser, M. E. Population growth in a wild bird is buffered against phenological mismatch. *Science* **340**, 488–491 (2013).

14. Møller, A. P., Rubolini, D. & Lehikoinen, E. Populations of migratory bird species that did not show a phenological response to climate change are declining. *Proc. Natl Acad. Sci. USA* **105**, 16195–16200 (2008).

15. Nylin, S. & Gotthard, K. Plasticity in life-history traits. *Annu. Rev. Entomol.* **43**, 63–83 (1998).

16. Roy, D. B. & Sparks, T. H. Phenology of British butterflies and climate change. *Glob. Chang. Biol.* **6**, 407–416 (2000).

17. Hodgson, J. A. et al. Predicting insect phenology across space and time. *Glob. Chang. Biol.* **17**, 1289–1300 (2011).

18. Roy, D. B. & Thomas, J. A. Seasonal variation in the niche, habitat availability and population fluctuations of a bivoltine thermophilous insect near its range margin. *Oecologia* **134**, 439–444 (2003).

19. Altermatt, F. Climatic warming increases voltinism in European butterflies and moths. *Proc. Biol. Sci.* **277**, 1281–1287 (2010).

20. Halupka, L., Dyrcz, A. & Borowiec, M. Climate change affects breeding of reed warblers *Acrocephalus scirpaceus*. *J. Avian Biol.* **39**, 95–100 (2008).

21. Van Dyck, H., Bonte, D., Puls, R., Gotthard, K. & Maes, D. The lost generation hypothesis: could climate change drive ectotherms into a developmental trap? *Oikos* **124**, 54–61 (2015).
22. Diamond, S. E., Frame, A. M., Martin, R. A. & Buckley, L. B. Species' traits predict phenological responses to climate change in butterflies. *Ecology* **92**, 1005–1012 (2011).

23. Thomson, L. J., Macfadyen, S. & Hoffmann, A. A. Predicting the effects of climate change on natural enemies of agricultural pests. *Biol. Control* **52**, 296–306 (2010).

24. Thackeray, S. J. et al. Trophic level asynchrony in rates of phenological change for marine, freshwater and terrestrial environments. *Glob. Chang. Biol.* **16**, 3304–3313 (2010).

25. Warren, M. S. et al. Rapid responses of British butterflies to opposing forces of climate and habitat change. *Nature* **414**, 65–69 (2001).

26. Garcia-Barros, E. Delayed ovarian maturation in the butterfly Hipparchia semeale as a possible response to summer drought. *Ecol. Entomol.* **13**, 391–398 (1988).

27. Stefanescu, C., Penuelas, J. & Filella, I. Effects of climatic change on the phenology of butterflies in the northwest Mediterranean Basin. *Glob. Chang. Biol.* **9**, 1494–1506 (2003).

28. Scali, V. Imaginal diapause and gonadal maturation of Maniola jurtica (Lepidoptera: Satyridae) from tuscany. *J. Anim. Ecol.* **40**, 467–472 (1971).

29. Visser, M. E. & Gienapp, P. Evolutionary and demographic consequences of phenological mismatches. *Nat. Ecol. Evol.* **3**, 879–885 (2019).

30. Suggitt, A. J. et al. Habitat associations of species show consistent but weak responses to climate. *Biol. Lett.* **8**, 590–593 (2012).

31. Roy, D. B. et al. Similarities in butterfly emergence dates among populations suggest local adaptation to climate. *Glob. Chang. Biol.* **21**, 3313–3322 (2015).

32. Lindestad, O., Wheat, C. W., Nylin, S. & Gotthard, K. Local adaptation of photoperiodic plasticity maintains life cycle variation within latitudes in a butterfly. *Ecology* **100**, e02550 (2019).

33. Gotthard, K., Nylin, S. & Wiklund, C. Adaptive variation in growth rate: life history costs and consequences in the speckled wood butterfly, Pararge aegeria. *Oecologia* **99**, 281–289 (1994).

34. Suggitt, A. J. et al. Extinction risk from climate change is reduced by microclimatic buffering. *Nat. Clim. Chang.* **8**, 713–717 (2018).

35. Oliver, T., Roy, D. B., Hill, J. K., Brereton, T. & Thomas, C. D. Heterogeneous landscapes promote population stability. *Ecol. Lett.* **13**, 473–484 (2010).

36. Pollard, E. A method for assessing changes in the abundance of butterflies. *Biol. Conserv.* **12**, 115–134 (1977).

37. Storkey, J. et al. Chapter One—The Unique Contribution of Rothamsted to Ecological Research at Large Temporal Scales, Vol 55. in *Advances in Ecological Research* (eds Dumbrell, A. J., Kordas, R. L. & Woodward, G.) 3–42 (Academic Press, 2016).

38. Fox, R., Warren, M. S., Brereton, T. M., Roy, D. B. & Robinson, A. A new red list of British butterflies. *Insect Conserv. Divers.* **4**, 159–172 (2011).

39. Heath, J., Pollard, E. & Thomas, J. A. *Atlas of Butterflies in Britain and Ireland* (Viking, 1984).

40. Fox, R. et al. Moths count: recording moths for conservation in the UK. *J. Insect Conserv.* **15**, 55–68 (2011).

41. Agassiz, D. J. L., Beavan, S. D. & Heckford, R. J. *Checklist of the Lepidoptera of the British Isles*. (Royal Entomological Society, 2013).

42. Menéndez, R., González-Megias, A., Lewis, O. T., Shaw, M. R. & Thomas, C. D. Escape from natural enemies during climate-driven range expansion: a case study. *Ecol. Entomol.* **33**, 413–421 (2008).

43. Farr, T. G. et al. The shuttle radar topography mission. *Rev. Geophys.* **45**, 1485 (2007).

44. Waring, P. & Townsend, M. *Field Guide to the Moths of Great Britain and Ireland* (Field Guides). (Bloomsbury Natural History, 2017).

45. Thomas, J. *Butterflies of Britain & Ireland*. (British Wildlife Publications, 2014).

46. U. K. Butterflies. Available at: http://www.ukbutterflies.co.uk/index.php. Accessed February 4, 2019.

47. UKMoths | Guide to the moths of Great Britain and Ireland. Available at: https://ukmoths.org.uk/. Accessed February 4, 2019.

48. Asher, J. et al. *Millennium Atlas of Butterflies in Britain and Ireland*. (OUP Oxford, 2001).

49. Brereton, T. M., Cruickshanks, K. L., Risely, K., Noble, D. G. & Roy, D. B. Developing and launching a wider countryside butterfly survey across the United Kingdom. *J. Insect. Conserv.* **15**, 279–290 (2011).

50. Moussus, J.-P., Julliard, R. & Jiguet, F. Featuring 10 phenological estimators using simulated data. *Methods Ecol. Evol.* **1**, 140–150 (2010).

51. Botham, M., Roy, D., Brereton, T., Middlebrook, I. & Randle, Z. United Kingdom Butterfly Monitoring Scheme: collated indices 2016. https://doi.org/10.5285/009334d6-ee33-43fc-a259-bf4d92fb69dc (2017).

52. Harrower, C.A., Bell, J.R., Blumgart, D., Botham, M.S., Fox, R., Isaac, N.J.B., Roy, D.B. & Shortall, C.R. UK moth trends from Rothamsted Insect Survey light-trap network (1968 to 2016). https://doi.org/10.5285/e7e0e4ad-f8c1-46fc-85b5-6d88057024b8 (NERC Environmental Information Data Centre, 2019).

53. Hollis, D. & McCarthy, M. UKCP09: Met Office gridded and regional land surface climate observation datasets (Centre for Environmental Data Analysis, 2017).

54. Kawahara, A. Y. & Breinholt, J. W. Phylogenomics provides strong evidence for relationships of butterflies and moths. *Proc. Biol. Sci.* **281**, 20140970 (2014).

55. R. Core Team. *R: A Language and Environment for Statistical Computing* (R Foundation for Statistical Computing, 2018).

56. Kearse, M. et al. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. *Bioinformatics* **28**, 1647–1649 (2012).

57. Wood, S. N. Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. *J. R. Stat. Soc. Ser. B Stat. Methodol.* **73**, 3–36 (2011).

58. Bates, D., Maechler, M., Bolker, B. & Walker, S. *lme4: Linear Mixed-effects Models Using Eigen and S4*. R package version 1.1-21 (2014).

59. Orme, D. et al. *Caper: Comparative Analyses of Phylogenetics and Evolution in R*. R package version 1.0.1 (2018).

60. Wickham, H. *ggplot2: Elegant Graphics for Data Analysis* (Springer, 2009).

61. Lucy, D. *blighty: United Kingdom Coastlines* R package version 3.1-4 (2012).

## Acknowledgements

We are grateful to all volunteers who contributed records to the four data sets used in this study. The UK Butterfly Monitoring Scheme is organised and funded by Butterfly Conservation, the Centre for Ecology and Hydrology, the British Trust for Ornithology and the Joint Nature Conservation Committee. The UKBMS is indebted to all volunteers who contribute data to the scheme. We are grateful to collaborators and staff who have contributed the data from the Rothamsted Insect Survey, a BBSRC-supported National Capability. Butterflies for the New Millennium and the National Moth Recording Scheme are run by Butterfly Conservation with funding from Natural England. We thank D. Boyes for assistance in classifying habitat specialisation for moths, E. Dennis and S. Freeman for helpful discussions surrounding the estimation of phenology using GAMs, K. Davis and A. Bakewell for advice on the construction of phylogenies, C. Shortall and A. Jansen van Rensburg for comments on the study design and results, and K. Das-mahapatra, G. Hurst and I. Owens for their involvement in designing the wider project. This work was supported by a grant from the Natural Environment Research Council (NERC; NE/N015797/1), and P.J.P. was supported by NERC grant NE/M013030/1.

## Author contributions

This study was instigated by C.J.M., C.D.T. and J.K.H. The study was primarily designed by C.J.M., C.D.T., D.B.R. and J.K.H. in discussion with M.A.B., J.R.Be., T.B., J.R.Br., R.F., G.M., I.M., P.J.P., R.R., I.S. and R.V. The statistical analysis was conducted by C.J.M., using the data provided by J.R.Be., R.F. and D.B.R.; and C.J.M. prepared the first draft of the paper. The study forms part of a wider programme of research which was originally designed by M.A.B., J.R.Be., J.R.Br., C.D., R.F., K.G., J.K.H., A.A.H., G.M., S.N., D.B.R., I.S., C.D.T. and C.W.W. All authors contributed substantially to revising the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

*Supplementary information* is available for this paper at https://doi.org/10.1038/s41467-019-12479-w.

*Correspondence and requests for materials should be addressed to C.J.M.*

*Peer review information* *Nature Communications* thanks Carol Boggs, Leslie Ries and Marcel Visser for their contribution to the peer review of this work. Peer reviewer reports are available.

*Reprints and permission information* is available at http://www.nature.com/reprints

*Publisher's note* Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_16.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution

4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2019