J A Hodgson

C D Thomas

T H Oliver

B J Anderson

T M Brereton

E E Crone

Predicting insect phenology across space and time Predicting insect phenology across space and time Predicting insect phenology across space and time Predicting insect phenology across

Keywords: climate change, emergence time, latitude, growing degree days, phenotypic plasticity, adaptation, butterfly monitoring scheme, lepidoptera climate change, emergence time, latitude, growing degree days, phenotypic 15 plasticity, adaptation, butterfly monitoring scheme, lepidoptera 16 17

Many species appear to be undergoing shifts in phenology, arising from climate change. To predict the direction and magnitude of future changes requires an understanding of how phenology depends on climatic variation. Species show large-scale spatial variation in phenology (affected by differentiation among populations) as well as variation in phenology from year-to-year at the same site (affected predominantly by local plasticity). Teasing apart spatial and temporal variation in phenology should allow improved predictions of phenology under climate change. This study is the first to quantify large-scale spatial and temporal variation in the entire emergence pattern of species, and to test the relationships found by predicting future data. We use data from up to 33 years of permanent transect records of butterflies in the UK to fit and test models for 15 butterfly species. We use generalised additive models (GAMs) to model spatial and temporal variation in the distribution of adult butterflies over the season, allowing us to capture changes in the timing of emergence peaks, relative sizes of peaks and/or number of peaks in a single analysis. We develop these models using data for 1973-2000, and then use them to predict phenologies from 2001-2006. For 6 of our study species, a model with only spatial variation in phenology is the best predictor of the future, implying that these species have limited plasticity. For the remaining 9 species, the best predictions come from a model with both spatial and temporal variation in phenology; for 4 of these, growing-degree-days have similar effects over space and

Global Change Biology

Many species appear to be undergoing shifts in phenology, arising from climate change. To 19 predict the direction and magnitude of future changes requires an understanding of how 20 phenology depends on climatic variation. Species show large-scale spatial variation in 21 phenology (affected by differentiation among populations) as well as variation in phenology 22 from year-to-year at the same site (affected predominantly by local plasticity). Teasing apart 23 spatial and temporal variation in phenology should allow improved predictions of phenology 24 under climate change. This study is the first to quantify large-scale spatial and temporal 25 variation in the entire emergence pattern of species, and to test the relationships found by 26 predicting future data. We use data from up to 33 years of permanent transect records of 27 butterflies in the UK to fit and test models for 15 butterfly species. We use generalised 28 additive models (GAMs) to model spatial and temporal variation in the distribution of adult 29 butterflies over the season, allowing us to capture changes in the timing of emergence peaks, 30 relative sizes of peaks and/or number of peaks in a single analysis. We develop these models 31 using data for 1973-2000, and then use them to predict phenologies from [2001][2002][2003][2004][2005][2006]. For 6 of 32 our study species, a model with only spatial variation in phenology is the best predictor of the 33 future, implying that these species have limited plasticity. For the remaining 9 species, the 34 best predictions come from a model with both spatial and temporal variation in phenology; for 35 4 of these, growing-degree-days have similar effects over space and time, implying high 36 levels of plasticity. The results show that statistical phenology models can be used to predict 37 phenology shifts in a second time period, suggesting that it should be feasible to project 38 phenologies under climate change scenarios, at least over modest time scales. 39 Many species appear to be undergoing shifts in phenology as a result of climate change 42 [START_REF] Sparks | The responses of species to climate over 2 centuries -an 514 analysis of the Marsham phenological record, 1736-1947[END_REF][START_REF] Crick | UK birds are laying eggs earlier[END_REF][START_REF] Menzel | Growing season extended in Europe[END_REF][START_REF] Roy | Phenology of British butterflies and climate change[END_REF]43 Stefanescu et al. 2003), and phenological changes constitute a high proportion of all of the 44 evidence that species are responding to climate change [START_REF] Walther | Ecological responses to recent climate change[END_REF]IPCC 2007). 45

Thanks to comprehensive biological recording, butterflies have become a model group in 46 detecting the impacts of climate change, including impacts on phenology (Sparks & Yates 47 1997;[START_REF] Roy | Phenology of British butterflies and climate change[END_REF][START_REF] Forister | Climatic trends and advancing spring flight of butterflies in 455 lowland California[END_REF][START_REF] Stefanescu | Effects of climatic change on the phenology of 519 butterflies in the northwest Mediterranean Basin[END_REF][START_REF] Dell | Climate change and the effect of increasing spring 447 temperatures on emergence dates of the butterfly Apatura iris (Lepidoptera : 448 Nymphalidae)[END_REF]48 Gordo & Sanz 2006;[START_REF] Menzel | Altered geographic and temporal 483 variability in phenology in response to climate change[END_REF]. However, understanding and predicting these 49 changes remains more of a challenge [START_REF] Visser | Keeping up with a warming world; assessing the rate of adaptation to 533 climate change[END_REF]. The climate change literature is replete 50 with statements about the potential detrimental effects of phenology shifts because they could 51 lead to significant changes in population growth rates and in community interactions (van 52 Asch et al. 2007;Doi et al. 2008). For example, temporal mismatches between plants and 53 their pollinators could reduce seed set in plants (including crops) and reduce food availability 54 for the pollinators [START_REF] Memmott | Global warming and the disruption of 477 plant-pollinator interactions[END_REF]). To predict whether such issues will become greater 55 in the future demands an ability to predict phenology shifts, based on climatic variation. 56

57

In insects, a variety of life history parameters (e.g. induction of diapause, size at pupation, 58 maternal investments, etc.) may evolve in order to maximize fitness under the prevailing 59 climate and the historic range of variability of that climate [START_REF] Davidson | The relationship between temperature and rate of development of insects at 445 constant temperatures[END_REF][START_REF] Lees | Diapause and photoperiodism in the fruit tree red spider mite 471 (Metatetranychus ulmi Koch)[END_REF]Manly 60 1974;[START_REF] Visser | Keeping up with a warming world; assessing the rate of adaptation to 533 climate change[END_REF]. Insects respond to very specific cues which, over their evolutionary 61 history, have helped them to stay in synchrony with important seasonal events, such as the 62 bud burst of trees [START_REF] Van Asch | Predicting adaptation of 528 phenology in response to climate change, an insect herbivore example[END_REF] [START_REF] Bale | Herbivory in global climate change research: 428 direct effects of rising temperature on insect herbivores[END_REF]. Because different species have evolved to use different 64 cues, and because climate change affects temperature but not day length, this could lead to 65 maladaptive responses (Doi et al. 2008). 66

67

A number of studies have used year-to-year variation in phenology (usually earliest 68 emergence or mean emergence date) regressed on yearly climate to predict future changes 69 (e.g. [START_REF] Roy | Phenology of British butterflies and climate change[END_REF][START_REF] Gordo | Temporal trends in phenology of the honey bee Apis mellifera (L.) 460 and the small white Pieris rapae (L.) in the Iberian Peninsula (1952-2004)[END_REF]. As well as temporal changes, many species 70 also show spatial variation in phenology (e.g. [START_REF] Zhou | Effects of 547 temperature on aphid phenology[END_REF][START_REF] Langvatn | Climate, plant phenology and 468 variation in age of first reproduction in a temperate herbivore[END_REF]White et 71 al. 1997), but a given climate variable may have different effects over space and over time 72 [START_REF] Pollard | Changes in the flight period of the hedge brown butterfly Pyronia-tithonus 496 during range expansion[END_REF][START_REF] Rock | Predicting geographical 505 and within-season variation in male flights of 4 fruit pests[END_REF][START_REF] Thompson | Spatio-temporal modelling and assessment of within-species 522 phenological variability using thermal time methods[END_REF]Doi et al. 2008;Doi & Takahashi 73 2008;[START_REF] Forkner | Timing is everything? Phenological 457 synchrony and population variability in leaf-chewing herbivores of Quercus[END_REF]. For example, the leafing date of japanese trees varies less with the 74 gradient of temperature over space than it does with changes in temperature over time at a 75 single site [START_REF] Doi | Latitudinal patterns in the phenological responses of leaf 452 colouring and leaf fall to climate change in Japan[END_REF]. If phenology varies to different extents in space and in 76 time this implies that species do not simply respond physiologically to temperature cues. 77

Instead they are likely to be, in part, locally adapted and/or responsive to spatially fixed cues, 78 such as daylength [START_REF] Bradshaw | Genetic shift in photoperiodic response correlated with 433 global warming[END_REF][START_REF] Visser | Keeping up with a warming world; assessing the rate of adaptation to 533 climate change[END_REF]. 79 80 Teasing apart the spatial variation in phenology, which is affected by local differentiation of 81 populations, and the temporal variation, which is affected by phenotypic plasticity within a 82 population, should allow much better predictions of phenology under climate change than 83 looking at one (space/time) dimension in isolation. However, few studies have considered 84 both spatial and temporal variation in phenology [START_REF] Rock | Predicting geographical 505 and within-season variation in male flights of 4 fruit pests[END_REF]Doi et al. 2008;Doi & 85 Takahashi 2008), and none of these has explicitly attempted to predict the future. 86 87

F o r R e v i e w O n l y 5
In this study we use data from up to 33 years of permanent transect records of butterflies in 88 the UK to fit and test models of the phenology of 15 butterfly species, modelling both spatial 89 and temporal variation in relation to the climate. We use generalised additive models (GAMs) 90 [START_REF] Wood | Generalized Additive Models: An Introduction with R. Chapman and Hall/ 543[END_REF] to model the distribution of adult butterflies over the season, independently of 91 the overall abundance at a particular site and year, allowing us to capture changes in the 92 timing of emergence peaks, relative sizes of peaks and/or number of peaks in a single 93 analysis. We develop these models using data for the time period 1973-2000, and

Butterfly data 100

We used butterfly abundance data from UK Butterfly Monitoring Scheme (collected and 101 maintained by Butterfly Conservation and the Centre for Ecology and Hydrology; Pollard & 102 Yates 1993;[START_REF] Asher | The Millenium Atlas of 426 Butterflies in Britain and Ireland[END_REF]. Butterflies of all species are recorded along a permanent 103 transect at weekly intervals between April and September at a network of sites across the 104 country. The number of sites taking part in the scheme has increased from 1 in 1973 to 750 in 105 2006. In order to fit models to roughly half the data, and use the other half to test model 106 predictions, we split the data set into 1973-2000 and 2001-2006 periods. 107 108 We focused on the 15 species that are known to have multiple generations per year in at least 109 part of their British range, since these have the most complex phenological patterns and might 110 be expected to be the most difficult to model and predict future changes. We ran all the 111 analyses separately for each species, using only sites where the species was found at a 112 reasonable frequency (usually using the condition that the median count over all weeks and 113 years should be greater than one, but this condition had to be relaxed for 4 species that were 114 particularly rare and/or temporally variable: holly blue, wall brown, small blue and wood 115 white; common English and Latin names are given in Table 2). 116 117

Climate data 118

We used historical monthly mean temperature data at the 10 km grid resolution, from 119 averages of the 5 km grid resolution data available from the Met Office ("UKCP09" data, 120 [START_REF] Met_Office | UKCP09: Gridded observation data sets 486[END_REF]. We decided to use growing degree days above 5° Celsius (henceforth 121 "GDD5") as a measure of the warmth available for development and reproduction that is 122

F o r R e v i e w O n l y
relevant to all butterfly species (even though there will also be various effects of winter 123 minima, summer maxima, rainfall, cloudiness, etc.) [START_REF] Romo | Use of degree-days in multiple-temperature experiments[END_REF]Bryant et al. 124 2002). From the monthly mean temperatures we calculated the GDD5 for each UKBMS site 125 (henceforth actual GDD5) in each year; defining "year" as the months between October of the 126 previous year (just after the previous transect season finished) and September (when the 127 transect season in question finished). Across the UK, GDD5 has been increasing in the time 128 period of this study, but there is also considerable year-to-year variation (Fig. 1). The average 129 increase in GDD5 between the period 1973-2000 and the period 2001-2006 is c.200, which is 130 nearly half of the inter-annual variation (variation in medians among years), and 10% of the 131 spatial variation in GDD5 that exists across all British 10km grid squares within a single year 132 (Fig. 1). 133 134Because we wanted to separate the effects of spatial variation in the climate from temporal 135 variation in the climate, we also calculated (a) mean GDD5 for each site over the period 1973-136 2000 (henceforth "site GDD5"), and (b) the difference between each site's actual GDD5 in 137 each year and the site mean (henceforth "GDD5 differential"). Hence each site and year has 3 138 interrelated GDD5 measures: (actual GDD5) t = (site GDD5) + (GDD5 differential) t where t 139 indexes years. 140 141

Statistical analyses 142

We used generalised additive models (GAMs) with a negative binomial error structure to 143 model the emergence pattern as a flexible, smooth function of the time of year (measured as 144 weeks from the first week in April). By using 2-D and 3-D smooths, we could represent the 145 fact that emergence pattern would change gradually with the climate over space and/or time. 146

We used the GAM implementation in the mgcv package in R [START_REF] Wood | Thin-plate regression splines[END_REF][START_REF] Wood | Fast stable direct fitting and smoothness selection for generalized additive 545 models[END_REF]; 147

F o r R e v i e w O n l y 8 R_Development_Core_Team 2009
) with the extension to negative binomial family in the 148 MASS package [START_REF] Venables | Modern Applied Statistics with S. Fourth Edition[END_REF]. Because the data consisted of butterfly counts (of 149 a certain species at a given site, week and year), but we wanted to investigate the emergence 150 pattern independently of the overall abundance of butterflies, we always included the average 151 abundance for the site and year as a linear covariate in the GAM. Specifically, this linear 152 covariate was calculated as: mean[ln(count i +1/27)] where i indexes the week, and there are 153 usually 26 weeks (sometimes fewer because of missed surveys). The mean is calculated in log 154 space because the negative-binomial GAM has a log link function, and a constant 1/27 is 155 added to the count because zeros cannot be logged. Using this average abundance as a 156 covariate means that the variation remaining to be explained is the distribution of individuals 157 among the weeks of a given year -which is the main quantity of interest -and that the 158 observation error structure relevant to such count data is retained. 159 160 We fitted 7 GAMs, representing different hypotheses about the drivers of phenology, to each 161 species' data for the early time period . We evaluated the GAMs on their ability to 162 predict the later time period (2001)(2002)(2003)(2004)(2005)(2006). Details of the 7 models used are given in Table 1. 163

The spatial component of climate is represented either by site GDD5 or by site northing: the 164 latter should be a better predictor of phenology if the species responds more to day length 165 The exact likelihood of each observation (x) given our model can be calculated as the 172 probability density at x of the negative binomial distribution with mean given by the GAM 173 prediction (according to the explanatory variables specific to that observation), and shape 174 parameter as fitted by the GAM (the same shape parameter for all observations under a given 175 model). Taking logs of these probabilities and summing over all observations (either 176 observations in the fitting period, or observations in the testing period) gave us the log 177 likelihood, L. The proportion of deviance "explained" was taken to be how far L lies between 178

Lnull (likelihood under model 00, see Table 1 were estimated for any model using the testing data). We can also use -2L to calculate Akaike 187 weights [START_REF] Burnham | Model Selection and Multimodel Inference: a Practical 438 Information-Theoretic Approach[END_REF] to give us a "degree of belief" in a certain model 188 compared to others. However, the Akaike weighting is strongly dependent on the assumption 189 that observations are independent -an assumption that was not really tenable for our data, 190 where there were multiple observations from exactly the same site. To show which of our 191 conclusions are robust to this assumption, we report raw model weights, and weights where 192 each observation only counts as 1/12 of a degree of freedom (L becomes L/12). The factor 12 193 was chosen based on the autocorrelation of residuals in the testing data. Residuals from 194 consecutive weeks at the same site and the same year had rank correlation coefficients >0.5, 195 but correlation dropped quickly with lag> 1 week, so roughly every other week's data could 196 that of 6 years of data, 5 may be pseudoreplicates. Dividing by 2 for non-independent weeks 199 and by 6 for non-independent years resulted in the factor 12. Although this "correction factor" 200 is crude, it should be noted that there is no perfect way of assessing independence, and that 201 the main results in this paper are robust to quite large changes in the weighting. 202 2, Fig. 2). For these 5 species, the model based on actual GDD5 (model 3) still 228 outperforms purely spatial models (models 1-2, Table 2, Fig. 2). This demonstrates that, even 229 though there may not be perfect space-for-time equivalence, there is a degree of consistency 230 in the response to temperature in different places (Table 2, Fig. 2). 231 232 For 9 species, northing is a consistently better spatial covariate than site average GDD5, i.e. 233 the Akaike weight of models with northing versus those with GDD5 is close to 1 (Table 2). 234

This suggests that day length could play a role in determining their phenologies, without 235 necessarily implying site-specific local adaptation. For 6 species, the opposite is true (low 236 weight of models with northing, compared to site average GDD5; Table 2); this suggests that 237 phenological changes are triggered by cues that are correlated with the average GDD5 of each 238 site, and this could potentially include population-specific photoperiod responses. Hypotheses 239 about population and species differences in photoperiod responses remain to be tested 240 experimentally. 241 242 There is great variation between species in exactly how phenology varies over space and time. 243

The number of emergence peaks appears to change in space and/or time for 7 out of the 15 244 species (Fig. 3, Fig. 4, Fig. S2). Six of these species show a greater number of peaks in the 245 southern (hotter) parts of their ranges than in the north (common blue (Fig. 3a), wood white, 246 small blue, small tortoiseshell, small copper, small heath ( S2), show higher late summer peaks in the south than in the north. The brown argus (Fig. 3d) 256

shows a higher first peak in the south. In the comma butterfly, the relative size of the middle 257 peak is larger in the south than in the north (Fig 4c-d). Finally, for some species the timing of 258 the peaks changes, with earlier emergences in the south and in hotter years (e.g. wall brown, 259 green-veined white, large white, Fig. S2). 260 Even though a large proportion of the deviance remained unexplained, predictions from our 264 models were surprisingly robust (similar amounts of deviance were explained in the testing as 265 in the training data). Because the data we fitted consisted of transect counts in particular 266 weeks and sites, which are inherently variable for many reasons, we would never expect 267 predictions to have a very high precision, but this should not prevent them from predicting 268 large-scale trends. The results show that statistical phenology models can be used to predict 269 phenology shifts in a second time period, suggesting that it should be feasible to project 270 phenologies under climate change scenarios, at least over modest time scales. 271 272 Our method isolates "phenology" in the sense of the distribution of adults over the season, 273 independently of the overall abundance in a particular site or year. This is a desirable first step 274 to understanding the underlying biology: it is more meaningful than other commonly used 275 metrics such as the date of first observation, mean date or date range, which can be biased by 276 changes in overall abundance or changes in voltinism. The entire pattern of emergence is 277 relevant to the interactions between species, such as between herbivores and plants or 278 One of the desirable features of our GAM modelling approach is its flexible ability to 284 represent the variety of phenological patterns seen in butterflies, especially multivoltine 285 butterflies. The patterns of emergence peaks in the GAMs (Fig. 3-4 andS2) are largely what 286 we would expect given what is known about the voltinism of these species. For example, the 287 common blue was already known to have univoltine populations in the north and bivoltine 288 ones in the south and this is picked up very clearly in the GAM (Fig. 3a). The northern peak is 289 timed midway between the two southern peaks, suggesting that it would not be simple for a 290 population to switch between the two modes over time, and thus is not surprising that the 291 northing-only model predicted best for this species. 292

293

The speckled wood also shows a complex pattern, broadly consistent with what is known 294 about its phenology. This species has a complex life cycle, overwintering in both larval and 295 pupal stages [START_REF] Nylin | Life-cycle regulation and life-history plasticity in 493 the speckled wood butterfly -are reaction norms predictable[END_REF][START_REF] Blakeley | Overwintering biology of Pararge aegeria[END_REF]). The pupae emerge somewhat before the 296 larvae in spring, producing two small peaks in the first half of the flight season (Fig. 3b). 297

Going from north to south, the two small peaks overlap more and more until they appear as 298 one extended peak. So we do not think this decrease in number of emergence peaks represents 299 a decrease in the numbers of generations the speckled wood achieves in a year -in fact it may 300 represent an increase. There could be 2 generations per year in the south but only 3 301 generations per 2 years in the north (the 3 rd peak being the offspring of the 1 st peak and the 2 nd 302 peak being the offspring of the 3 rd peak in the previous year). It is a little surprising that the 303 northing-only model gave the best predictions for this species, because it has been observed 304 elsewhere that development times are very dependent on temperature, and that the flight 305 period can be extended (start earlier and finish later) in relatively warm years at a single site 306 [START_REF] Shreeve | The effect of weather on the life cycle of the speckled wood butterfly 512 Pararge aegeria[END_REF]). This highlights the potential for more detailed species-specific models to 307 peaks in the south than in the north. For the holly blue and the green-veined white, the second 313 peak in the north is smaller than the first, and may represent only a partial second generation. 314

For the large white the second peak is always substantially bigger than the first, so it seems 315 more reasonable to assume that there are 2 complete generations per year, but that high 316 breeding success and/or survival during summer are higher in the warm south. The brown 317 argus (Fig. 3d) shows a higher first peak in the south, seemingly because the first adult peak is 318 extremely small in northern locations rather than because summer breeding success is 319 increased in the north (although the latter is possible because of reduced parasitism in the 320 north; [START_REF] Menendez | Escape from 479 natural enemies during climate-driven range expansion: a case study[END_REF]). In the comma butterfly, the relative size of the middle peak is 321 larger in the south than in the north (Fig 4c-d). Two types of comma adults are produced in 322 the mid-summer peak: bright orange, active individuals that breed again in the same year, 323 with their offspring emerging in late summer; and darker, less active individuals that will 324 overwinter directly as adults and that will not breed until the following year [START_REF] Nylin | Effects of changing photoperiods in the life-cycle regulation of the Comma 488 butterfly, Polygonia-c-album (Nymphalidae)[END_REF]325 1992). Therefore, the shift most likely reflects an increase in the frequency of two 326 generations per year in the south. 327 328

Insight into the drivers of phenology 329

Our statistical models can potentially capture the net effects of phenological adaptations (the 330 average and slope of the reaction norms) that differ between populations, and make successful 331 shown for many insects, that phenological events are triggered by a combination of 333 temperature and day-length cues. For example, the developmental pathways of the comma 334 butterfly, for which model 5 (indicating an interaction effect of latitude and GDD5 335 differential) was the best predictor, are known to be affected by both photoperiod and 336 temperature [START_REF] Nylin | Effects of changing photoperiods in the life-cycle regulation of the Comma 488 butterfly, Polygonia-c-album (Nymphalidae)[END_REF][START_REF] Nylin | Seasonal plasticity in life-history traits -growth and development in 490 Polygonia-c-album (Lepidoptera, Nymphalidae)[END_REF]. Any ongoing evolution of the reaction norms will cause 337 mismatches between our models and future data. To predict this evolution is a very difficult 338 task, and would require more detailed knowledge of the traits that underlie the realised 339 phenology [START_REF] Visser | Keeping up with a warming world; assessing the rate of adaptation to 533 climate change[END_REF]. However, our results suggest testable predictions about which 340 species respond to which triggers. 341

342

The relative performance of our 5 alternative models suggest that different species currently 343 have different amounts of regional adaptation of phenology, and different levels of plasticity 344 in phenology. Some species show quite flexible space-for-time substitution (i.e., the 345 phenology can simply be predicted by knowing the temperature accumulated in a given 346 site/year), but others show patterns of emergence that are less flexible over time. Six out of 347 the 15 species' phenologies were best fit solely by fixed spatial effects (latitude or average 348 GDD5 of a site), four species showed flexible space-for-time substitution, while five species 349 showed complex patterns that were best fitted by an interaction between spatial and temporal 350 variation in phenology (i.e. change in phenology over time was dependent on the spatial 351 location of the site). Species with fixed spatial effects have, prima facie, no plasticity and 352 therefore may only respond to climate change gradually, through evolutionary changes (for 353 example evolution of critical photoperiods that induce diapause, Bradshaw & Holzapfel 354 2001). When species exhibit space-for-time equivalence this suggests that reaction norms are 355 the same in all populations, and that no evolution is required to respond to moderate climate 356 change. Therefore we may expect these species to "keep up" with climate change, and we 357

F o r R e v i e w O n l y
18 might also have greater confidence in the predictions of our models for future years for these 358 species. When species exhibit space-time interactions, this suggests that reaction norms are 359 different in different geographic areas. Therefore we may expect these species to show an 360 intermediate rate of response to climate change, or to respond in some geographic areas but 361 not others. It is not clear whether year-to-year plasticity necessarily increases fitness, 362 compared to a fixed phenology strategy -presumably this would also depend on the 363 phenological responses of the resources, mutualists and natural enemies of the species under 364 consideration. 365 366 Limitations 367

We have focused on how the GDD5 summed over a whole year affects the entire sequence of 368 emergence, and in doing so we have inevitably glossed over many of the climatic and 369 biological details that underlie these phenomena. Better fits could undoubtedly be obtained in 370 future studies by including additional climate variables, particularly those that show little 371 correlation with GDD5. The small tortoiseshell in particular could be an interesting subject 372 for further study because so much variation was left unexplained by temperature or northing. 373

The speckled wood deserves further consideration because its phenology responds to early 374 and late season temperatures [START_REF] Shreeve | The effect of weather on the life cycle of the speckled wood butterfly 512 Pararge aegeria[END_REF]) that were not fully encapsulated within our 375 GDD5 approach, and the same may be true for additional species. It may be argued that better 376 predictions of the future would be obtained from a mechanistic model of the species' 377 development, life-history and population dynamics. In several experimental studies it has 378 been shown that degree-days have quite a direct relationship to the development rate of each 379 larval instar of insects [START_REF] Manel | Modeling insect development time of two or more larval stages 473 in the field under variable temperatures[END_REF][START_REF] Gu | Statistical estimation of degree days of mosquito development 463 under fluctuating temperatures in the field[END_REF]Thompson & Clark 380 2006;[START_REF] Raworth | Initiation of oviposition after winter diapause in the spider mite 502 Tetranychus urticae (Acari : Tetranychidae): prediction and historical patterns[END_REF][START_REF] Trnka | European 525 Corn Borer life stage model: Regional estimates of pest development and spatial 526 distribution under present and future climate[END_REF]). These relationships are often used to build 381 sophisticated models of the phenology of certain well-known species, such as crop pests (e.g. 382 [START_REF] Collier | Accumulated temperatures for predicting the time of emergence 440 in the spring of the cabbage root fly, Delia radicum (L) (diptera, Anthomyiidae)[END_REF][START_REF] Trnka | European 525 Corn Borer life stage model: Regional estimates of pest development and spatial 526 distribution under present and future climate[END_REF]. However, such models usually need to be calibrated 383 differently for different geographic areas, reflecting the variability of development conditions 384 in the field relative to laboratory conditions, microclimate differences among sites and 385 individuals, and the genetic differentiation of populations [START_REF] Collier | Accumulated temperatures for predicting the time of emergence 440 in the spring of the cabbage root fly, Delia radicum (L) (diptera, Anthomyiidae)[END_REF]Bryant et al. 386 2002;[START_REF] Trnka | European 525 Corn Borer life stage model: Regional estimates of pest development and spatial 526 distribution under present and future climate[END_REF]. Because physiological data will be required for large numbers of 387 populations of every species, it is unlikely to be feasible to parameterize such models for the 388 majority of species for which we would like to predict phenological responses. 389 390 Whilst our model projections performed well in testing data, the predictions were only for a 391 relatively short testing period (2001)(2002)(2003)(2004)(2005)(2006), immediately following the model training period 392 , and the climate during the testing period remained largely within the bounds of 393 conditions experienced during the training period (Fig. 1). As such, the testing data were 394 independent, but we remain cautious about the capacity of such models (or other models) to 395 project phenologies far into the future. Reaction norms may evolve, especially once the new 396 climate falls outside the range of the conditions previously experienced by those populations, 397 and cause future phenologies to deviate systematically from model projections. This would 398 be true even for species that already have high levels of phenotypic plasticity in phenology. 399

F o r R e v i e w O n l y 19
For multivoltine species, such as these butterflies, extra generations may be observed under 400 future, hotter climates, which could accelerate the pace of demographic and evolutionary 401 responses. 402 403 Although we remain cautious about long-term projections, the approach we have taken 404 appears suitable to model and predict short-term phenology patterns in space and time. The 405 statistical models presented here can be parameterized from the sort of phenology data that 406 are regularly collected by volunteer recorders, and this approach can straightforwardly be 407 the overall mean for the period of fitting data and the period of testing data 558 (2001)(2002)(2003)(2004)(2005)(2006). 559 compared to GAM 0 (where phenology is a smooth function of week of the year, with no 561 spatial or temporal variation). Bars are shaded to aid the comparison of similar models: light 562 grey when site GDD5 is the spatial covariate, dark grey when northing is the spatial covariate, 563 and white when actual GDD5 is the covariate (combining both spatial and temporal 564 variation). For the absolute amounts of deviance explained, see supplementary material, Fig. 565 S1. 566 GAM included temporal as well as spatial variation. In each panel, the unlogged, GAM-fitted 577 distribution of butterflies across the 26 weeks of transect surveys is shown by lines, 578 standardized so that the sum of all weeks equals 1. Different lines represent different values of 579 Northing: equally-spaced steps between the southernmost (red) and the northernmost (blue) 580 parts of the species' data (the range, in British National Grid10 km squares is also printed in 581 the panel). To illustrate temporal variation, predictions for an average year (GDD5 582 differential= 0) and a hot year (GDD5 differential= 200), are shown in adjacent panels. For 583 example, the small copper exhibits two distinct peaks in mid-late summer in the hot south, but 584 only in hot years. Equivalent figures for all other species are given in supplementary material, 585

  high levels of plasticity. The results show that statistical phenology models can be used to predict phenology shifts in a second time period, suggesting that it should be feasible to project phenologies under climate change scenarios, at least over modest time scales

  than temperature. The temporal component is represented by GDD5 differential. If the basic 166 responses to temperature are the same everywhere (i.e. different populations are not adapted 167 to the climate differently) then the actual GDD5 experienced by each site in each year might 168 be the best predictor of phenology.

  ) and Lsat. Lsat is the maximum likelihood 179 possible for those data: that which would occur under a model where predictions are exactly 180 the same as observations (i.e. a saturated model) and where a Poisson distribution (no 181 overdispersion) is assumed. So proportion deviance explained = 1-(L-Lsat)/(Lnull-Lsat). 182 183 Our results focus on comparing the log likelihood of the testing (2001-2006) data under each 184 GAM. Using a separate test data set effectively penalizes overfitted models, and the model 185 with the highest log likelihood can be said to be the best model (i.e., no GAM parameters 186

  pseudoreplicate. Residuals from the same week at the same site for different 197 years had rank correlation coefficients >0.5, irrespective of the time lag in years -suggesting 198

  variation in GDD5 to determine observed emergence patterns (model 4 or 227 5 best; Table

  Fig S2)). Two of these species 247 show an increased tendency to produce an extra peak in the south in relatively hot years 248 (small copper (Fig. 4a-b), small heath (Fig S2)). The speckled wood (Fig 3b) appears to be an 249 exception in that there are fewer emergence peaks in the far south (2) than further north (3). 250 of the complex life cycle of this species, fewer peaks might not imply 251 fewer generations per year (see discussion). 252 253 For other species, the number of peaks does not seem to change, but the relative sizes of the 254 different peaks do change. The holly blue (Fig 3c), green veined white, and large white (Fig 255

  first study to quantify large-scale spatial and temporal variation in the entire 262 emergence pattern of species, and to test the relationships found by predicting future data. 263

  between plants and pollinators. Observing the distribution of adults over the season in relation 279 to habitat and climate factors also permits fuller consideration of how climate change alters 280 the interactions among species, and how it might change selection pressures.

  , the GAMs pick up changes in the relative sizes of the different peaks. The 311 holly blue (Fig 3c), green veined white, and large white (Fig S2), show higher late summer 312

  predictions without necessarily knowing what the climate triggers are. It is likely, as has been 332

  numbers of species. The models have a predictive power which is reasonably 408 high and surprisingly consistent across species. We see this as an important step towards 409 anticipating which species might be threatened as a result of phenological mismatches as the

Figure 1 :

 1 Figure 1: the distribution of GDD5 across all 10 km squares in the UK (boxplots of the 556

Figure 2 :

 2 Figure 2: Comparison of deviance in an independent test dataset explained by GAMs 1-5 560

Figure 3 :Figure 4 :

 34 Figure 3: Spatio-temporal variation in phenology of selected species where the best-predicting 567

Fig

  Fig. S2. 586

  Figure2

  . Temperature and day length (photoperiod) may both 63

	be used as cues
	F o r
	R
	e v
	i e
	w
	O n l
	y
	4

y Methods 99

	94	then use
	95 96 97	each model to predict phenologies for the period 2001-2006. The testing data enables us to assess our ability to predict detailed phenological patterns into the future. Furthermore, the relative predictive power of alternative models that have different predictor variables allows F o r
	98	us to make inferences about the factors that regulate phenology in different species. R
		e v
		i e
		w
		O n l

Tables 551 Table 1 :

 5511 summary of 7 GAM models fitted to each species' data 552

	Code number	Model variables †	Max. # of	Interpretation
			parameters
	00	abundance	2	No climate effects and no seasonality
	0	00+ smooth(week)	12	No climate effects
	1	00+ smooth(week * site gdd5)	25	Spatial variation only (suggests response to fixed
				site-specific cues that are correlated with site
				temperature; possibly local photoperiod
				adaptation)
	2	00+ smooth(week * northing)	25	Spatial variation only (suggests response to fixed
				cues that are correlated with latitude; probably
				photoperiod)
	3	00+ smooth(week * actual gdd5)	25	Space-for-time substitution (suggests plastic
				response to temperature)
	4	00+ smooth(week * site gdd5 * gdd5 diff )	41	Different spatial & temporal effects (suggests site-
				specific responses, as in model 1, combined with
				flexible temperature responses)
	5	00+ smooth(week * northing * gdd5 diff)	41	Different spatial & temporal effects (suggests
				likely photoperiod response, as in model 2,
				combined with flexible temperature response)
	† each site and year has 3 interrelated GDD5 measures: (actual GDD5) t = (site GDD5) + (GDD5 differential) t . Site GDD5 is the mean
	GDD5 for the site for the period 1973-2000.

Table 2 :

 2 summary of results for each species 553

species common English name Latin name no. obervations in 1973-2000 period no. obervations in 2001-2006 period best % explained in fitting data (1973-2000) best % explained in testing data (2001-2006) best model at predicting testing data second best model at predicting testing data Akaike weight of model 3 relative to 1 & 2* Akaike weight of 3 with df correction* Akaike weight of the models with northing vs. those with site GDD5* Akaike weight of models with northing with df correction* Akaike weight of model 0 relative to all others* Akaike weight of 0 with df correction*

	Large white	Pieris brassicae	20281 15059	20	20	3	4	1.00000 1.00000 0.04052 0.43445 0.00000 0.00000
	Adonis blue	Polyommatus	1030	591	36	27	3	1	0.99997 0.60770 0.00000 0.18239 0.05358 0.05358
		bellargus
	Wall brown	Lasiommata megera	35070 15463	25	20	3	5	1.00000 0.99991 1.00000 0.99950 0.00000 0.00000
	Holly blue Common blue Speckled wood Pararge aegeria Celastrina argiolus Polyommatus icarus Green-veined white Pieris napi Small tortoiseshell Aglais urticae Small heath Coenonympha pamphilus Small copper Lycaena phlaeas Comma Polygonia c-album F o r R 38157 33980 21900 17347 42782 33303 35484 20578 14272 8252 F o r R e v i e w 18 15 1 23 23 2 18 16 2 20 19 4 10 4 4 20111 12892 22 21 5 7184 3087 22 18 5 2021 2538 21 18 5 *model weights are by definition between 0 and 1, but are only shown to 5 decimal places 0 5 3 5 3 4 4 3 O n l y 0.00000 0.00003 0.00000 0.00260 0.08465 0.08465 0.00000 0.00000 1.00000 1.00000 0.00000 0.00000 0.00000 0.00299 1.00000 0.99998 0.00000 0.00000 1.00000 1.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.99984 0.00000 0.00000 0.00000 0.00000 0.99988 0.64557 1.00000 0.93065 0.00000 0.00000 1.00000 0.90559 1.00000 0.91856 0.00001 0.00001 1.00000 0.77540 1.00000 0.75212 0.00335 0.00335 e v
	Brown argus Small blue Wood white	Aricia agestis Cupido minimus Leptidea sinapis	2488 877 i 2151 933 e 911 803 w	24 37 37	22 30 23	2 2 2	0 5 3	0.00000 0.00011 1.00000 0.88283 0.12504 0.12504 0.00001 0.18354 0.99999 0.76765 0.06272 0.06272 0.08686 0.31741 0.99784 0.53666 0.07144 0.07144
	Small white	Pieris rapae	29420 20736 O 20 n l 19	3	4	1.00000 1.00000 0.00000 0.25855 0.00000 0.00000 30
					y

Acknowledgements