This is the accepted version of the article: Carnicer i Cols, Jofre; Stefanescu, Constantí; Vives Ingla, Maria; [et al.]. Phe- notypic biomarkers of climatic impacts on declining insect populations : A key role for decadal drought, thermal buffering and amplification effects and host plant dynamics. DOI 10.1111/1365-2656.12933 This version is avaible at https://ddd.uab.cat/record/203312 under the terms of the license 1 Phenotypic biomarkers of climatic impacts on declining insect populations: a key role for 2 decadal drought, thermal buffering and amplification effects and host plant dynamics 3 4 Jofre Carnicer1,2*, Constantí Stefanescu2,3, Maria Vives-Ingla1, Carlos López2, Sofia Cortizas1, 5 Christopher Wheat4, Roger Vila5, Joan Llusià2, Josep Peñuelas2 6 7 1 Department of Evolutionary Biology, Ecology and Environmental Sciences, University of Barcelona, 8 08028 Barcelona, Spain 9 2 CREAF, Global Ecology Unit, Autonomous University of Barcelona, 08193 Cerdanyola del Vallès, 10 Spain 11 3 Natural History Museum of Granollers, E-08402 Granollers, Spain 12 4 Department of Zoology (Population Genetics), University of Stockholm, Svante Arrheniusväg 13 18B, Stockholm University, S-10691 Stockholm, Sweden 14 5 Institute of Evolutionary Biology (CSIC-UPF), 08003 Barcelona, Spain 15 16 Corresponding author: jofre.carnicer@ub.edu 17 18 Running head: Phenotypic biomarkers of climate impacts and insect population declines. 19 20 Keywords: climate change, butterflies, phenotypic biomarker, multiannual drought, thermal 21 buffering, host plant, Pieris napi. 22 23 Type of paper: primary research article. 24 25 26 27 1 1 Abstract 2 3 1. Widespread population declines have been reported for diverse Mediterranean 4 butterflies over the last three decades, and have been significantly associated to 5 increased global change impacts. The specific landscape and climatic drivers of these 6 declines remain uncertain for most declining species. 7 2. Here we analyse whether plastic phenotypic traits of a model butterfly species (Pieris 8 napi) perform as reliable biomarkers of vulnerability to extreme temperature impacts in 9 natural populations, showing contrasting trends in thermally exposed and thermally 10 buffered populations. 11 3. We also examine whether improved descriptions of thermal exposure of insect 12 populations can be achieved by combining multiple information sources (i.e. integrating 13 measurements of habitat thermal buffering, habitat thermal amplification, host plant 14 transpiration, and experimental assessments of thermal death time (TDT), thermal 15 avoidance behaviour (TAB) and thermally induced trait plasticity). These integrative 16 analyses are conducted in two demographically declining and two non-declining 17 populations of P. napi. 18 4. The results show that plastic phenotypic traits (butterfly body mass and wing size) are 19 reliable biomarkers of population vulnerability to extreme thermal conditions. Butterfly 20 wing size is strongly reduced only in thermally exposed populations during summer 21 drought periods. Lab rearing of these populations documented reduced wing size due to 22 significant negative effects of increased temperatures affecting larval growth. We 23 conclude that these thermal biomarkers are indicative of the population vulnerability to 24 increasing global warming impacts, showing contrasting trends in thermally exposed 25 and buffered populations. 26 5. Thermal effects in host plant microsites significantly differ between populations, with 27 stressful thermal conditions only effectively ameliorated in mid-elevation populations. 2 1 In lowland populations we observe a six-fold reduction in vegetation thermal buffering 2 effects, and larval growth occurs in these populations at significantly higher 3 temperatures. Lowland populations show reduced host plant quality (C/N ratio), 4 reduced leaf transpiration rates and complete aboveground plant senescence during the 5 peak of summer drought. Amplified host plant temperatures are observed in open 6 microsites, reaching thermal thresholds that can affect larval survival. 7 6. Overall, our results suggest that butterfly population vulnerability to long-term drought 8 periods is associated to multiple co-occurring and interrelated ecological factors, 9 including limited vegetation thermal buffering effects at lowland sites, significant 10 drought impacts on host plant transpiration and amplified leaf surface temperature, as 11 well as reduced leaf quality linked to the seasonal advance of plant phenology. Our 12 results also identify multi-annual summer droughts affecting larval growing periods as a 13 key driver of the recently reported butterfly population declines in the Mediterranean 14 biome. 15 16 17 18 19 20 3 1 Introduction 2 Declines in butterfly populations across diverse species over the last three decades have been 3 described in the Mediterranean basin (Stefanescu et al., 2004, 2011ab; Wilson et al., 2005, 4 2007; Carnicer et al., 2012, 2013; Zografou et al., 2014; Melero et al., 2016). Negative effects 5 of land use changes and global warming have been proposed as the main drivers of the observed 6 declining trends (Stefanescu et al., 2004, 2011ab; Wilson et al., 2005, 2007). These negative 7 demographic trends affect both habitat generalist and specialist butterfly species in the 8 Mediterranean biome, and the spatial diversity of most functional groups (e.g. host-plant use, 9 dispersal capacity, habitat specialisation and thermal niche) is negatively associated with 10 increased temperatures and aridity (Stefanescu et al., 2011ab, Carnicer et al., 2013). 11 Furthermore, the available evidence suggests that global-warming induced population responses 12 are intimately linked to complex interactions with habitat features and host plant dynamics 13 (Merrill et al., 2008; Suggitt et al., 2012; De Frenne et al., 2013; Oliver et al., 2014, 2015; 14 Carnicer et al., 2017). In line with this idea, it has been suggested that specific habitat attributes 15 can modify global warming impacts on butterfly populations, triggering both positive and 16 negative demographic responses. For example, it has been shown that the densification of forest 17 habitats associated with land abandonment can cool local microclimates, buffering the impacts 18 of global warming in some plant and insect populations and resulting in positive or neutral 19 demographic responses to global warming (De Frenne et al., 2013; Nieto-Sánchez et al., 2015). 20 On the other hand, populations inhabiting sites lacking effective habitat thermal buffering could 21 experience increased negative impacts of extreme temperatures, resulting in substantial long- 22 term demographic declines (Parmesan, et al. 2000). In addition to the effects of habitat thermal 23 buffering, the thermal exposure of butterfly populations can be crucially determined by other 24 key processes, such as the seasonal variation of host plant transpiration and leaf water content 25 during summer drought, the operation of thermal amplification processes in microhabitats or the 26 display of thermal avoidance behaviours in the insect larvae allowing the selection of cool 27 microsites at the host plant (Carnicer et al. 2017). These key co-acting processes are often not 4 1 measured and their complex interactions remain poorly described. To understand the relative 2 importance of all these processes, integrative studies combining multiple information sources in 3 intensively studied populations are warranted. 4 5 Here we provide an integrative study of the thermal exposure in four populations of Pieris napi, 6 combining multiple sources of information (demographic and climatic data, phenotypic trait 7 data, measurements of habitat thermal buffering, host plant traits, and experimental assessments 8 of thermal responses). Furthermore, we explore whether temperature-responsive phenotypic 9 traits can be applied as reliable biomarkers of the different vulnerability to increased 10 temperatures in these intensively studied populations. Ample experimental evidence supports 11 that diverse life history and functional traits of butterflies are highly responsive to temperature 12 variation and show predictable responses to extreme temperature treatments (Jones et al., 1982; 13 Sheridan & Bickford, 2011; Bauerfeind & Fischer, 2013ab, 2014; Nail et al., 2015). In 14 particular, wing and body size measures have been identified as traits highly responsive to 15 temperature variation and climate change impacts (Atkinson, 1994; Atkinson & Sibly, 1997; 16 Nygren et al., 2008; Kingsolver, 2009; Talloen et al.; 2009; Sheridan & Bickford, 2011; Forster 17 et al., 2012). Therefore, it is likely that an extensive quantification of plastic phenotypic traits in 18 declining and non-declining natural populations could indicate their different vulnerability to 19 warmer conditions. In other words, if a specific morphological trait of a species is known to 20 respond plastically and in a linear manner to thermal conditions, then we can potentially deduce, 21 for specific populations, the exposure to these thermal conditions by quantifying its 22 morphology. Moreover, if we measure extreme thermal conditions in a target population, which 23 should induce a negative morphological response, and find non-altered biomarker values, we 24 can suspect that the population is buffered from stressful conditions by microhabitat effects. 25 In the Mediterranean region, summer drought periods and increased summer temperatures are 26 tightly linked and significantly associated (Fig. S1). Therefore, during extreme summer drought 27 periods we expect phenotypic traits to be affected by extreme temperature impacts. In this 5 1 context, those population sites lacking effective habitat thermal buffering effects should present 2 a significant negative response in temperature-responsive biomarker traits. In contrast, we 3 expect that populations characterised by effective microsite buffering mechanisms should 4 present non-significant trends in temperature-responsive phenotypic biomarkers (see 5 supplementary text S1 for a formal definition of the term biomarker and a simple mathematical 6 framework supporting this definition). 7 To test this hypothesis and develop an integrative analysis of thermal exposure in a butterfly 8 species, we address the following five research objectives using the green-veined white Pieris 9 napi as a species model: i) to analyse whether plastic phenotypic traits perform as reliable 10 biomarkers of vulnerability to extreme temperature impacts in natural populations, by 11 comparing phenotypic trait responses in four populations of Pieris napi; ii) to experimentally 12 estimate thermal death time responses (thermal susceptibility (z) and critical thermal limit 13 (CTmax)) and the thermal threshold for avoidance behaviour (TAB) for this model species, iii) to 14 quantify thermal buffering in microsites, assessing whether they provide non-stressful thermal 15 habitats only in specific localities; iv) to evaluate whether host plant resource dynamics 16 qualitatively differ between the studied populations, and v) to assess whether increased drought 17 impacts could explain the reported long-term population declines in the selected model species. 18 19 Methods 20 Study species 21 Pieris napi is a widely distributed Holarctic butterfly, common across most of North America 22 and Europe, though only locally in North Africa. Throughout its distribution, it shows a clear 23 preference for humid habitats, such as wetlands, riparian forests and irrigated agricultural land. 24 In Catalan lowland areas, there is a succession of 4-5 generations from early spring (March- 25 April) to autumn (October-early November), with overwintering in the pupal stage. Maximum 26 abundance is typically recorded in early summer, in coincidence with the peak of the third 6 1 generation. This peak is followed by a period of 1-2 months when abundance is much reduced, 2 in coincidence with summer drought. Butterflies then reappear by the end of September, in what 3 normally constitutes the last annual generation. In mountain areas, where the phenology is 4 constrained by colder temperatures, a succession of three generations from April to September 5 is the most common pattern. At most montane sites, abundance increases all over the season and 6 reaches its maximum in the third and last annual generation. Eggs are laid singly on a wide 7 range of wild Brassicaceae, Lepidium draba and Brassica nigra being the two most common 8 host plants in lowland areas, and Alliaria petiolata, Arabis glabra and Cardamine pratensis 9 those mostly used in mountain habitats. Other secondary host plants have been recorded over 10 the region (García-Barros et al., 2013). 11 Study zone 12 We studied two lowland declining populations (sites 1 and 2) and two mid elevation non- 13 declining populations (sites 3 and 4) in Catalonia, NE Spain. Sites 1 and 2 were located at two 14 protected coastal wetlands (Delta del Llobregat and Aiguamolls de l’Empordà, 133 km apart). 15 In contrast, sites 3 and 4 were located at higher elevations, also in natural protected areas (Zona 16 Volcànica de la Garrotxa (503 m a.s.l.) and Montseny (1031 m a.s.l.), 41 km apart). Mid 17 elevation sites were characterised by a heterogeneous mosaic of different habitat types, 18 including open fields, small wetland and riverine areas, and temperate and evergreen forests. A 19 more detailed summary of the geographic, climatic and ecological attributes of the selected sites 20 is provided in Table S1 and Fig. S2. To quantify long term demographic trends, sites were 21 surveyed from 1994 to 2012 as part of the Catalan Butterfly Monitoring Scheme 22 (www.catalanbms.org) via weekly butterfly counts along fixed transect routes from March to 23 September (a total of 30 recording weeks per year). All individuals seen within 2.5 m on each 24 side and 5m in front of the recorder were counted, using the standard methodology of the 25 Butterfly Monitoring Schemes (Pollard & Yates, 1994; Schmucki et al., 2016). For site 1 26 demographic surveys were available only for seven years distributed in two discrete periods 7 1 [1994-1997, 2007-2009]. An annual index calculated as the sum of weekly counts was used as 2 the measure of population abundance at each season. 3 4 Population trends, climate and landcover data and model selection approach 5 Climate factors and landscape use changes have been identified as the main drivers of long-term 6 butterfly population trends in the Mediterranean biome (Stefanescu et al., 2011ab). However, 7 detailed models combining climatic and dynamic landscape data are still warranted to 8 quantitatively assess the relative contribution of these two factors to long-term butterfly 9 demographic declines. For this purpose, we compiled a database integrating butterfly annual 10 abundance indices, monthly climatic rainfall and temperature data for the 1994-2012 time 11 period (Domingo-Marimón, 2015), and landcover dynamics data for 1994-2012. To study 12 landscape dynamics, aerial orthoimages (1:25000) for 1993, 2001, 2006 and 2012 were 13 digitised using MiraMon, a geographic information system (Pons, 2002). The images were 14 provided by the Cartographic Institute of Catalonia (http://www.icgc.cat/en/). We selected a 15 circular area (2 km of diameter) around the field transect sites and quantified the changes in the 16 total surface (m2) of the following nine landcover types: wetland and continental water (L1), 17 dense forest (L2), sparse forest (L3), shrubland (L4), grassland and herbaceous meadows (L5), 18 urbanised land (L6), bare land (L7), road/lane areas (L8) and beach area (L9). A continuous 19 annual sequence of estimated land cover changes for 1994-2012 was obtained applying spline 20 fits using JMP (SAS Institute, 2012) and saving predicted values between consecutive 21 orthoimages in the time series. 22 To analyse the observed temporal trends in the butterfly annual abundance of the four 23 populations over 1994-2012 spline fits and ordinary least squares models (OLS) were 24 implemented. Two model selection approaches were sequentially applied, first using only 25 climatic variables (approach 1) and subsequently combining land and climatic variables in an 26 integrated model (approach 2). The first modelling approach was simply used to reduce the 27 large number of climatic variables analysed (a total of 28 monthly temperature and precipitation 8 1 variables). In other words, we first selected monthly climatic variables significantly associated 2 with the observed butterfly demographic trends (OLS step-wise approach, (SAS Institute, 3 2012)) and then we combined the selected climatic variables and dynamic landscape data in an 4 integrated model selection approach. All possible models computable in each approach were 5 contrasted in terms of their corrected Akaike’s Information Criterion (AICc) and Bayesian 6 Information Criterion (BIC), and the models with the lowest values were selected. The 7 explanatory power of competing variables was contrasted by the stepwise selective approach 8 and by comparing the estimates for the selected predictors (JMP package, SAS Institute, 2012). 9 Digitised orthoimages were not available for site 4, precluding the inclusion of this site in the 10 landscape modelling analyses. We included in the model selection approach monthly climatic 11 variables of two consecutive years in order to account for the previous autumn growing period 12 of winter diapausing generations (i.e. current and previous year climatic data). 13 14 Phenotypic biomarkers of population vulnerability 15 In order to identify phenotypic traits that could perform as climate-extreme biomarkers, 16 butterfly populations were intensively sampled with weekly resolution during 2014 and 2015, 17 covering the whole flying period (early spring to late autumn). Weekly samples were composed 18 of a minimum of four males and four females. Supplementary samples were collected during 19 seasonal abundance peaks. A total of 1265 butterflies were finally collected (see Table S2 for 20 further details). The following phenotypic traits and their seasonal variance were quantified in 21 the four selected populations: dry body mass, dry wing mass, wing size (i.e. length and area). 22 Dry wing mass and wing area variables were significantly correlated and were considered 23 synonymous descriptors (R2=0.40; p<0.0001). The same was true for dry body mass and wing 24 area measures (R2=0.41; p<0.0001). In addition, we also quantified wing melanism, and whole 25 body ∂13C, ∂15N, %N and %C. However, these variables were not strongly related to climate 26 variability and were discarded. 9 1 To quantify wing size (length and area), wing samples were photographed using standardised 2 settings (fixed Nikon D7100 with a SigmaMacro objective at a height of 41.5 cm). The 3 quantification was performed with ButterflyPhotoGUI, a Matlab algorithm (developed by 4 Hedrick, T.; Kingsolver lab, University of North Carolina), so that wing size corresponded to 5 the number of pixels in the area defined by three fixed landmarks in the hindwing (tip of the 6 vein M1, tip of the vein CuA1, and the intersection of the veins CuA1-CuA2, Fig. S3). The wing- 7 vein naming system applied is described in Wahlberg et al. (2014). In addition, wing length 8 (mm) and area (mm2) were measured independently in a subset of standardised photographs. All 9 these measures (wing area in pixels, wing area in mm2 and wing length in mm) were strongly 10 correlated and, thus, were considered related descriptors (R2>0.80; p<0.0001). Wing size was 11 finally chosen to present our graphical results (i.e. wing length 1 in Fig. S3). 12 To track the impacts of extreme summer climatic conditions on wing size, we focused on the 13 weekly variation of this trait during the spring, summer and autumn periods. To more precisely 14 quantify the effect of climatic variables on phenotypic trait variability, we modelled the 15 variation of wing size using ordinary least squares models (OLS) and introducing the following 16 predictor variables: site, year (2014/2015), mean temperature during the larval and pupal growth 17 period (25 days previous to the adult collection), accumulated rainfall previous to adult 18 collection (60 days), mean relative humidity previous to adult collection (60 days), photoperiod 19 (mean of 25 days previous to adult collection), and sex (male/female). The interactions between 20 the predictor variables were examined and significant interactions were retained. For each 21 climatic variable, different possible temporal spans were assessed, ranging from 5 to 120 days 22 (with a 5-day resolution), compiling subsets of related climatic variables. The climatic variables 23 that were finally selected were characterised by higher correlation coefficients with the 24 modelled variables (wing size) in multivariate correlation analyses (JMP, SAS Institute, 2012). 25 We excluded the first generation, i.e. winter-diapausing individuals, from the modelling 26 analyses. Photoperiod was calculated following Kirk (1994). We randomly collected both 27 freshly emerged and older, worn individuals, and estimated adult age by quantifying wing 10 1 condition state using an ordinal scale (Fig. S5). No significant effects of wing condition were 2 observed on wing size models 3 4 Common-garden experiment 5 To assess whether the selected biomarker traits were reliable predictors of direct temperature 6 effects on the phenotype, we performed a common garden split-brood experiment. Five female 7 lines from the lowland site 2 were initiated, with eggs reared under photoperiod conditions 8 inducing direct development (13:11 L:D). The offspring (eggs) were divided in two temperature 9 treatments (20 ºC, 25ºC). As illustrated in Fig. S6, the 20 ºC treatment corresponds to the 10 observed mean daily temperatures in June or late August. In contrast, the 25 ºC treatment 11 corresponds to the warmest mean daily temperatures of June, July and August in the study 12 period (Fig. S6, and see Bauerfeind & Fischer 2013ab, 2014 for additional experimental studies 13 in this species). 14 Fresh leaves of the host plant species Lepidium draba were provided to the larvae ad libitum. A 15 total of 143 adult individuals belonging to 5 different families were finally obtained (see Table 16 S3 for detailed numbers). The experiment allowed testing the effects of treatment, sex and 17 family on wing size. The heritability of the measured traits was estimated using MCMCglmm 18 (Hadfield 2010; Aalberg Haugen et al., 2012). To assess putative differences between 19 populations in plastic phenotypic responses between low-elevation and mid-elevation sites, 32 20 additional adult individuals, belong from mid-elevation site 3 were assessed in replicated 21 experimental split-brood conditions (20/25 ºC treatments, two female lines). 22 23 Host plant microsite climatic measures 24 Maximum summer temperatures often surpass the critical thermal limits of invertebrate 25 ectotherms in multiple biomes (Sunday et al., 2014), and as a result a key role of thermal 26 buffering processes has been identified for population persistence (e.g. Ashton et al., 2009; 11 1 Sunday et al., 2014; Suggitt et al., 2015; Pateman et al., 2016). Consequently, a robust 2 evaluation of climate impacts on butterfly populations requires quantifying microclimatic 3 thermal variability and habitat buffering effects at the host plant level during larval growth 4 periods. In addition, the analysis of the temporal variation of host plant traits over the season 5 allows the identification of critical periods of resource scarcity and changes in host plant 6 quality. In the studied system, the dominant host plants were Lepidium draba at lowland 7 wetland areas and Alliaria petiolata at mid elevation sites. For L. draba and A. petiolata six host 8 plant microsites were selected at lowland and mid elevation mountain ranges, respectively (see 9 Table S4 for details). In each microsite we installed an automatic temperature and humidity 10 sensor (LascarElectronics EL-USB-2-LCD) recording hourly climatic variability over 2014 and 11 2015. In the lowlands, L. draba was mostly distributed in open microsites (open meadows and 12 grassland areas) and more rarely under tree canopy and/or shrub cover. Egg-laying by P. napi 13 has been recorded on plants growing in both conditions. Four microsite sensors were therefore 14 distributed in the most representative open meadow microsites, and two sensors were located at 15 closed-canopy host plant microsites to quantify the effect of canopy cover on temperature and 16 humidity records. An additional and commonly used host plant, Brassica nigra, was also 17 present in lower numbers at lowland site 2 (inhabiting open microsites, along ditches). Two 18 sensors were located at Brassica nigra microsites to quantify the observed trends for this host 19 plant. At mid elevation sites, a single dominant host plant (A. petiolata) was preferentially 20 located and used for egg-laying at closed-canopy sites. However, a comparatively smaller 21 number of host plants were distributed in open meadow and/or grassland microsites. Four 22 sensors were located at the dominant and representative conditions (closed-canopy sites) and 23 two additional sensors were located at the more unusual open grassland microsites. To contrast 24 microsite climatic measures and standard measures, daily temperature and rainfall records were 25 obtained from four meteorological stations located nearby the four surveyed transects (2-5 km) 26 and at the same altitudinal range (Table S5). The automatic temperature and humidity sensors 27 (LascarElectronics EL-USB-2-LCD) were located at 25 cm height above the soil surface using 12 1 metal stakes, and were protected from direct solar radiation by a plastic envelope sustained by a 2 wire-mesh cylinder (installed 5 cm above the sensor and thus precluding the direct incidence of 3 solar radiation). The sensors were surrounded by the host plant leaves and were also covered by 4 abundant herbaceous vegetation (the herbaceous layer ranged between 50-120 cm of height). At 5 closed sites, the sensors were in addition directly affected by the shadows of the surrounding 6 shrubs and trees. The sensors estimated the air temperature, relative humidity and dew point 7 with hourly resolution. 8 9 Thermal avoidance behaviour (TAB) and thermal death time (TDT) experiments 10 To avoid the exposure to critical thermal temperatures, butterfly larvae may display 11 thermal avoidance behaviors (i.e. short movements to cooler microsites of the host 12 plant). However, the thermal thresholds for these behaviours remain poorly quantified 13 and experimentally studied in most butterfly species. We experimentally assessed 14 thermal avoidance behaviour in 149 larvae of Pieris napi (last instar, site 2), assessing a 15 total thermal range of 28-48 ºC. Larvae were firstly placed on a leaf of a potted Alliaria 16 petiolata, and acclimated at an ambient temperature of 22 ºC for 5 minutes. Then, we 17 experimentally raised leaf surface temperatures to a selected Celsius degree treatment 18 (in the thermal range of 28 - 48 ºC) using a 70 W light lamp and carefully controlling 19 the leaf surface temperature with a HANNA HI935005N thermal sensor. Larval 20 behaviour was recorded for 2 minutes, annotating three types of responses: thermally 21 neutral, thermally positive and thermal avoidance movements to cooler microsites. For 22 each Celsius degree treatment, we assessed 5-7 larvae. Each larvae was used in a single 23 thermal trial. A two parameter logistic model was fitted to model the changes in the 24 frequency of thermal avoidance behaviour (f(T)) with increasing leaf surface 25 temperature (T): 13 1 1 f (T) = 1+ e-aT-b 2 where a is the growth rate of the function and b is the thermal inflection point (ºC) in 3 which we observed a 0.5 frequency of thermal avoidance behaviours. 4 Thermal death time experiments (TDT) allow predicting from first principles when 5 environmental temperatures may affect larval survival (Deutsch et al 2008, Rezende et 6 al 2014). To assess the upper critical thermal limit (CTmax) in Pieris napi we 7 implemented a static thermal death time experiment with three static thermal treatments 8 (40, 42 and 44 ºC; Rezende et al. 2014). We estimated CTmax and thermal susceptibility 9 (z) from the equation: 10 Tko = CTmax - zlog10 t 11 where t is the observed time to death of last instar Pieris napi larvae in static thermal 12 experimental treatments, Tko corresponds to the constant stressful temperature levels 13 applied, CTmax is the temperature that would result in knowckdown or death at 1 min 14 (log10 t = 0) and z is the constant of thermal susceptibility describing how thermal 15 tolerance decays with the duration of the heat challenge. The experiment was 16 implemented in 60 individuals from 6 family lines collected at site 2. 20 individuals 17 were assessed in each thermal treatment (40, 42 and 44 ºC). 18 19 Plant trait measurements 20 To evaluate whether host plant resource dynamics qualitatively differed between populations we 21 quantified the weekly variation of leaf %N and leaf C/N. Previous empirical works have shown 22 that nitrogen strongly determines butterfly host plant quality of mature leaves (Mattson, 1980; 23 Slansky & Feeny, 1977; Scriber & Slansky, 1981; Myers, 1985; Kaitaniemi et al., 1998). 24 Moreover, leaves containing less nitrogen constrain insect performance and reduce pupal mass 14 1 in field experiments (Myers, 1985; Kause et al., 1999, but see Fischer & Fiedler 2000). Plant 2 phenology and drought have been identified as key drivers of leaf nitrogen variation (Kause et 3 al., 1999; Grant et al., 2014). Notably, drought and phenology should produce qualitatively 4 different effects on the selected nitrogen-related traits (%N, C/N). In the case of phenology, a 5 progressive reduction on the quantity of nitrogen in the leaves should be expected with plant 6 maturation and leaf ontogeny due to the translocation of nitrogen-rich resources to flowers, 7 fruits and rhizomes (Kause et al., 1999; Jacobs, 2007). Overall, phenology should promote a 8 progressive reduction of leaf %N and an increase of leaf C/N over the late spring and early 9 summer period. In contrast, in the case of drought-induced effects on plant ecophysiology, an 10 increase in the quantity of nitrogen in the leaves could be expected (e.g. Bauerfeind & Fischer 11 2013b; Grant et al., 2014; Valim et al., 2016). To perform the weekly leaf measurements five 12 plants were sampled per week and site, collecting 5-8 leaves per plant for the analyses (Table 13 S6). Seasonal trends for C/N ratio and %N were modelled applying ordinary least squares for 14 linear trends and spline fits for non-linear trends (JMP package (SAS Institute 2012)). 15 In order to achieve a more detailed assessment of plant responses to drought stress at lowland 16 sites, monthly measures of leaf stomatal conductance and leaf surface temperatures were 17 specifically conducted for plants L. draba and B. nigra (using a LICOR 6400 portable 18 photosynthesis system). These measurements were restricted to open microsites and were 19 conducted at midday (12.00-14.00) (Aiguamolls de L’Empordà wetlands, site 2). Four 20 replicates were measured for each host plant species. The measurements were conducted in 21 2015, a year characterised by warm and dry summer conditions (Fig. S7). Stomatal conductance 22 and leaf surface temperatures have been widely applied as integrated indicators of drought and 23 heat stress in herbaceous plants (Munns et al., 2010; Anissa et al., 2013). Experimental 24 evidence shows that leaf conductance and surface temperature show qualitatively different 25 responses in heat treatments and soil drought experiments (Anissa et al., 2013). In conditions of 26 abundant soil moisture, Brassica plants respond to strong air temperature stress with an increase 27 in leaf stomatal conductance values, producing in turn positive leaf-to-air temperature 15 1 differences (i.e. cooler leaf temperatures relative to air temperatures due to increased leaf 2 transpiration (Anissa et al., 2013)). In contrast, if plants experience combined soil water and air 3 heat stress, which is probably more common in Mediterranean ecosystems, reduced leaf 4 conductance values and negative leaf-to-air temperature differences are to be expected (i.e. 5 higher leaf temperatures relative to air temperatures due to reduced leaf transpiration). 6 Photosynthetic rates (A) and stomatal conductances (gs) were measured between 12:00 a.m. and 7 14:00 p.m. at a quantum flux density (PPFD) of 1080±19 µmol m-2 s-1 and ambient air 8 temperature under a controlled CO2 concentration of 400±2 ppm. To conduct the measurements 9 one leaf was enclosed in a clamp-on gas-exchange cuvette of 2 cm2. We selected healthy leaves 10 that were not affected by insect larvae consumption and/or fungal damages. Air flow through 11 the dynamic cuvette was 732±0.05 ml min-1. A Licor-6400XT (4647 Superior Street P.O. Box 12 4425 Lincoln, Nebraska USA) gas-exchange system was used. 13 14 Results 15 Climatic and population trends 16 Model selection analyses using only climatic data identified June rainfall of the current year as 17 the best predictor of the interannual variation of Pieris napi abundance (Tables S7 and S8). An 18 analysis of the temporal trends for this climatic predictor (June rainfall) over 1994-2012 19 identified a decadal period of increasing drought (1997-2007, Fig. 1a). The observed decadal 20 reduction in June rainfall was highly significant at the four sites (Fig. 1a). In line with the 21 reported decadal trend of increasing summer drought stress, butterfly abundance at lowland sites 22 significantly declined, paralleling the trend of June rainfall (Fig. 1b). As a result, lowland 23 populations showed a sharp reduction of more than one order of magnitude respect the initial 24 abundance numbers. In contrast, mid elevation site populations remained fairly stable over the 25 1997-2007 drought period (Fig 1c), and were therefore not paralleling June rainfall trends as 26 observed at the lowland sites. After the decadal drought period, however, population at site 3 27 increased significantly (R2=0.38, p=0.0024, Fig. 1c), and this increase was significantly 16 1 correlated to an abrupt increase in June rainfall during 2008-2012 (R2=0.38, p=0.0028). 2 Lowland sites showed significantly lower June rainfall values during the 1997-2007 period 3 (Tukey-Kramer test, R2=0.35; p=0.0005, Fig. S8). 4 5 6 Figure 1. Climatic and butterfly demographic dynamics over the 1994-2012 period. a) Annual 7 variation of June rainfall at the four population sites. Significant linear rainfall trends are 8 indicated in the 1997-2007 period. b) Observed variation of butterfly annual abundance at 9 lowland sites (Delta del Llobregat (1) and Aiguamolls de l’Empordà (2) protected wetlands). A 10 significant polynomial and a spline fit are illustrated. c) Observed variation of butterfly annual 11 abundance at mid elevation sites (Zona Volcànica de la Garrotxa (3) and Montseny Ranges (4) 12 protected areas). Spline and linear fits were applied. When significant, the variance explained 13 by the linear fit (R2) is indicated. 17 1 2 The model selection approach combining climatic and landcover data reported that both types of 3 variables significantly contributed to the reported demographic trends. Overall, however, the 4 estimates of the models suggested a stronger and predominant effect for climatic variables in the 5 reported trends (June rainfall, Tables S7–S11). For landscape variables, significant negative 6 effects of reduced meadow cover extent during 1994-2012 were detected at site 2. At mid 7 elevation site 3, a positive effect of increased wetland area was detected. Landscape data were 8 not available for 2013-2015 and therefore these years were excluded from the butterfly annual 9 abundance models. Nevertheless, we examined the observed population trends for an extended 18 1 period (1994-2015) and the results were fully consistent with the trends reported for 1994-2012 2 (Fig. S9). Mid elevation site 4 showed a stable population trend (Fig. 1 C), and no significant 3 relationships with climatic variables in OLS models were observed for this site. Population 4 abundances at low elevation sites areas were higher (Fig. 1), presumably due to a much higher 5 spatial density of host plants per unit of surface observed in these wetland areas. 6 7 Different population sensitivity to temperature impacts 8 The analysis of the environmental variation of wing size revealed significantly different trends 9 between lowland and mid elevation populations (Fig. 2a). Of note, these trends were 10 significantly associated to increased summer temperature only at lowland sites (Fig. 2a). 11 Consistent with this, OLS models identified the seasonal variation of temperature during the 12 larval and pupal growth period (25 days previous to the adult emergence and collection) as the 13 principal driver of diverging wing size seasonal trends and detected site x temperature 14 interactions (Table 1). The interactions between site and temperature were highly significant, 15 reporting strong negative effects only at lowland sites (Table 1, and see Fig. 2a). 16 Consistent with the field observations, the split brood common-garden experiment demonstrated 17 a significant link between temperature and wing size variation (Fig 2b). Tables S12-S15 18 summarise the results of the split-brood experiment. Significantly different wing size values 19 were observed for the 20 ºC and the 25 ºC treatment, with reduced wing lengths observed for 20 the high temperature treatment. In addition, we observed significant effects of sex and family, 21 with females showing significantly lower wing sizes (Tables S12-S15). The effect of the 22 temperature treatment, however, was dominant and stronger than family and sex effects. Wing 23 size heritability estimate reported by MCMCglmm models was 0.41 (CI=0.12-0.76). 24 Importantly, the observed wing sizes for the 25 ºC treatment were consistent with the observed 25 range of wing size values in the field dataset in the same range of temperatures (25 ºC [i.e. mean 26 temperature during the growth period (25 days)], Fig. 2a, grey square area). Significant effects 27 of 20/25 ºC thermal treatment on larval developmental times were observed, as reported in Fig. 19 1 S10 (R2=0.71, p<0.0001). Site effects were not significant (Tables S14-S15). Mid altitude sites 2 showed more plastic responses to temperature in experimental treatments (Table S15), 3 indicating that the flat trends in Fig 2a were not related to a lack of thermally induced wing 4 plasticity in mid altitude populations. 5 6 Figure 2. Observational (a) and experimental (b) trends in the selected phenotypic biomarker 7 trait (wing size (mm)). a) Observed relationships between butterfly wing size and environmental 8 temperature at the four sites. The grey square represents the mean wing size and the 95% 9 confidence interval of the thermal stress treatment (25 ºC), matching the field observational 10 values (lines) at lowland sites. b) Observed differences in wing size measurements between two 11 experimental temperature treatments (20 / 25 ºC). The line across each diamond represents the 12 treatment mean. Diamond plots indicate the 95% confidence interval for each treatment (vertical 13 span) and mean (midpoint line). Green dots represent mid-elevation individuals. Red dots 14 represent low-elevation individuals. We concluded that experimental and field results were in 15 agreement, suggesting a key role of stressful temperatures at lowland sites in the reported wing 16 size trends. c) Estimated thermal inflection point for behavioural avoidance responses in last 17 instars of Pieris napi. d) Observed thermal death time (TDT) in static thermal treatments in last 18 instar larvae of Pieris napi. The line within the box represents the median sample value. The 19 ends of the box represent the 25th and 75th quantiles. 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 21 1 2 3 Table 1. OLS model of the variation of wing size. Values in bold highlight the principal effect 4 of temperature variation and site*temperature interactions (negative in lowland, declining 5 populations; and positive at mid elevation sites). 6 Wing size Model fit: R2=0.38, p<0.0001, AICc=18232.3, BIC=18300.8 Estimate Std Err t p Intercept 317368.16 190999.1 1.66 0.0970 Temperature -9164.388 1321.999 -6.93 <.0001 (Temp) Site 1 -13424.82 8506.542 -1.58 0.1150 Site 2 -23883.71 8205.615 -2.91 0.0037 Site 3 32545.087 10085.57 3.23 0.0013 Site 4 4763.4443 6388.543 0.75 0.4561 Sex (female) -14702.02 2308.956 -6.37 <.0001 Year -153.3687 4089.377 -0.04 0.9701 Photoperiod 10464.797 5677.805 1.84 0.0657 Rainfall 764.4488 7598.182 0.10 0.9199 Relative Humidity 2546.9183 2090.97 1.22 0.2236 Temp*site1 -6824.945 1318.059 -5.18 <.0001 Temp*site2 -7078.373 1352.495 -5.23 <.0001 Temp*site3 4663.1566 1517.146 3.07 0.0022 Temp*site4 9240.1615 1799.712 5.13 <.0001 Temp*sex -1274.876 854.2187 -1.49 0.136 (female) 7 8 9 22 1 Thermal avoidance behaviour (TAB) and thermal death time (TDT) experiments 2 The results of the thermal avoidance behaviour experiment are shown in Figure 2c. The 3 observed behavioural response of thermal avoidance was well described by a two parameter 4 logistic model (R2=0.92, p<0.0001) with the following parameters: growth rate a = 0.86±0.24, 5 thermal inflection point b =37.46±0.37 ºC. These results indicate a rapid shift to behavioural 6 avoidance responses at temperatures above 37.5 ºC in the last instar larvae of Pieris napi. 7 Fig. 2d synthesises the results of the static thermal death time experiments. Thermal death time 8 experiments (TDT) reported an estimate of the temperature resulting in death at 1 min of 9 exposure (CTmax) of 51 ºC, and a thermal susceptibility constant (z) of 4.11±0.33 (ºC). The 10 observed TDT relationships for 100% and 50% of mortality are illustrated in Fig. S11a. The 11 thermal threshold for a time of exposure equal to the whole larval period was estimated in 32.5 12 ºC (Fig S11b). For a daily exposure of 6 hours to maximum daily temperatures (TE6h, 10 am- 13 16 pm), the TDT curve indicates a thermal threshold of 34.5 ºC (Fig. S11b). These thermal 14 thresholds were achieved in warm summer days characterised by mean daily temperature ≥ 25 15 ºC in 2014-2015 (Fig. S11c). 16 17 Thermal stress during summer drought and microsite effects 18 Analysis of meteorological data during 2014-2015 for the four population sites found that mean 19 June maximum temperatures were in the range of 28-30 ºC at three sites (1-3) and around 25 ºC 20 at site 4 (Fig. 3a). Maximum daily temperatures reached the experimental TDT thermal 21 thresholds in warm summer days (i.e. for thermal values higher than the quantile 75th, Fig 3a). 22 Next, we examined whether host plant microhabitat buffering effects at the four sites could 23 allow reduced maximum temperature values. Analysis of host plant microsite climatic data 24 (2014-2015) revealed strong buffering effects only at mid elevation populations (-5.2±0.17 ºC), 25 and only for those plants located at closed microsites (Fig. 3b, green rectangles; Table S16). In 26 contrast, in low elevation sites we observed limited cooling effects of canopy cover at closed 23 1 sites (-0.79±0.32 ºC, Fig. 3b, Table S16). Open microsites were characterised by amplified 2 mean maximum June temperatures (1.9±0.18 ºC in lowland sites and 3.04±0.28 ºC in mid 3 altitude sites; orange rectangles in Fig. 3b). The analysis of daily cycles of temperature variation 4 showed that temperatures of warm summer days reached values higher than the experimental 5 TDT thresholds for several hours in open sites (Fig. 3c). In contrast, these range of thermal 6 values were not achieved in closed microsites of mid elevation sites (most values < 30 ºC, Fig. 7 3d). Overall, we conclude that significantly different thermal buffering effects were observed at 8 lowland and mid elevation sites, in line with the previously reported trends for butterfly 9 demographic declines and for phenotypic biomarkers (wing size responses). 10 11 Host plant resource dynamics 12 The analysis of the seasonal patterns of host plant availability and quality (C/N ratio) revealed 13 important differences between mid elevation and lowland sites. Mid elevation sites were 14 characterised by a continuous availability of fresh Alliaria petiolata leaves during the whole 15 summer period and, consequently, by more stable temporal C/N ratios (Fig. 4a). In contrast, at 16 lowland sites, the leaves of the two host plants Lepidium draba and Brassica nigra presented a 17 significant linear increase in the C/N ratio (indicating a progressive reduction of host plant 18 quality with the advance of summer and plant phenology). This trend culminated in total leaf 19 senescence at the end of June – early July (Fig. 4b). 20 21 22 23 24 25 26 27 24 1 Figure 3. Comparison of June maximum temperature measurements (i.e. mean of the daily 2 maximum temperatures during June) at standardised meteorological stations and at host plant 3 microsites for 2014-2015. a) Meteorological data. The line within the box represents the 4 median. The ends of the box represent the 25th and 75th quantiles. The lines that extend from 5 the box indicate the following distances: 25th quantile - 1.5*(interquartile range) and 75th 6 quantile + 1.5*(interquartile range). The plane yellow dotted line indicates a thermal threshold 7 of 32.5 ºC calculated from the TDT relationship, corresponding to a time of thermal exposure 8 equivalent to whole larval period (LP, see Fig. S11). The orange dotted line indicates a 34.5 ºC 9 threshold, corresponding to the TDT for 6 hours of daily exposure to maximum temperatures 10 over the larval period (TE6h). The red dotted line indicates the experimental threshold for 11 thermal avoidance behaviour (TAB) of 37.5 ºC. Different capital letters indicate significantly 12 different means (Tukey-Kramer test). The grey surface area illustrates the logarithmic decrease 13 of the thermal death time with linearly increasing temperatures. b) Temperature-humidity host 14 plant sensor data. Green rectangles indicate the observed habitat buffering effect in Celsius 15 degrees at host plant microsites relative to standardized meteorological records. Orange 16 rectangles indicate the observed thermal amplification of host plant microsites relative to 17 standardized meteorological records. c and d) Observed daily variation of June temperatures at 18 two host plant microsites characterised by contrasting buffering trends (c, open microsite, 19 lowland site; d, closed microsite, mid elevation site). A spline fit (black line) indicates the mean 20 trend observed. A smooth surface illustrating the density of data points is provided. Red contour 21 lines indicate maximum point density. The contour lines are quantile contours in 5% intervals 22 (i.e. 5% of the temperature measurements are below the lowest (blue) contour, 10% are below 23 the next contour. The highest (red) contour has about 95% of the thermal values below it). 24 25 1 2 3 Therefore, in contrast to mid elevation sites, summer drought at the lowland sites produced a 4 relatively large period (45-65 days) in which we observed a total absence of fresh leaves due to 5 the complete senescence of the aboveground organs (leaves and shoots, corresponding to Julian 6 days 190 (July the 9th) - 235 (August the 23rd) in Fig. 4b). In agreement with these observations, 7 a significant reduction of leaf conductance values was observed during summer drought (Fig. 8 5a). The observed values were below 0.2 mol/m2 s (Brassica nigra: 0.176±0.034 mol/m2 s; 9 Lepidium draba: 0.137±0.021 mol/m2 s). These values matched the range of conductance 10 values reported in water stress experiments for Brassica species in stressful conditions (Anissa 11 et al., 2013; Guo et al., 2015). Complementary results for photosynthetic rates (A) and sub- 12 stomatal CO2 concentrations (ci) for Lepidium draba are reported in Fig. S12. 13 In accordance with these trends, significantly higher temperatures at the leaf surface in relation 14 to air temperature were recorded during the peak of summer drought (Tukey-Kramer test, L. 15 draba, R2=0.64, p<0.0001; B. nigra, R2=0.56, p<0.0001, Fig. 4c and d). Similarly, with the 26 1 onset of summer season, midday leaf temperatures significantly increased (L. draba temperature 2 increase: TJune – TMay = 14.5 ºC; B. nigra, TJune – TMay = 14.4 ºC; Tukey-Kramer test, p<0.0001). 3 As a result, midday leaf temperatures in June reached 37.56±0.35 ºC in Brassica nigra and 4 38.16±0.52 in Lepidium draba (Fig. 5c and d). In summary, during the peak of summer drought 5 the results for lowland plants indicated significant reductions on leaf quality (increased C/N 6 ratios), significantly reduced conductance values (gs< 0.2 mol/m2 s) and significantly increased 7 leaf surface temperatures. 8 9 Figure 4. Observed annual variation of leaf host plant quality (C/N content ratio). Higher C/N 10 ratio corresponds to lower host plant quality. a) Observed trends at mid elevation site 3. b) Low 11 elevation sites 1 and 2. The grey surface area illustrates the summer drought period. 12 13 14 15 27 1 2 Figure 5. Observed monthly variation of leaf conductance for lowland plants Lepidium draba 3 (a) and Brassica nigra (b). The line across each diamond represents the mean. The vertical span 4 of each diamond represents the 95% confidence interval. Leaf conductance is linked to plant 5 transpiration and leaf energy balance, and hence to the ability of the plant to cool itself under 6 heat stress. c and d) Observed midday June temperatures of the leaf surface of low elevation 7 host plants Lepidium draba (upper panel) and Brassica nigra (lower panel). Temperatures were 8 significantly higher than air temperatures synchronously recorded using LICOR 6400 portable 9 photosynthesis system (indicated as “Air” in the panels). Expected thermal death time for Pieris 10 napi is provided in the right axis for the amplified leaf surfaces temperatures. The line within 11 the box represents the median sample value. The ends of the box represent the 25th and 75th 12 quantiles. The grey surface area illustrates the logarithmic decrease of the thermal death time 13 with linearly increasing temperatures. 14 15 16 17 18 28 1 2 Discussion 3 Our results indicate that a decadal trend of increased summer drought has triggered long-term 4 declines of P. napi populations at lowland sites (Fig. 1). In contrast, the analysis of 5 microclimatic conditions experienced by mid elevation populations suggests a key role for 6 habitat buffering processes in these sites. Mid elevation populations presented significantly 7 stronger thermal buffering effects in closed habitat microsites (a six fold increase), and as a 8 result individuals were characterised by comparatively larger butterfly wing sizes than 9 butterflies from lowland populations exposed to similar conditions of environmental thermal 10 stress without access to thermal buffering micro-refugia (Figs 2 & 3). Moreover, a continuous 11 availability of high quality leaf resources with low C/N ratios was observed only at mid 12 elevation sites, in non-declining populations (Fig. 4). In contrast, lowland host plants showed a 13 progressive seasonal reduction of leaf nitrogen content, possibly associated to the seasonal 14 advance of the flowering and fruiting phenological cycle (Fig. S13). In line with this finding, 15 previous studies have documented a decrease in crude protein and digestible fibre after 16 flowering in Lepidium draba (Jacobs, 2007). Finally, consistent effects of summer drought were 17 observed in the leaf conductance of Lepidium draba and Brassica nigra, the two lowland host 18 plants, resulting in turn in significantly increased leaf surface temperature. In addition to these 19 combined drought and leaf heat stress impacts, lowland populations also experience periods of 20 host plant resource scarcity (during late July-August), which are in turn associated with a 21 decrease in abundance due to pupal aestivation (Fig. 4). Overall, our results suggest that 22 butterfly population vulnerability to long-term drought periods is associated to multiple co- 23 occurring and interrelated ecological factors, including limited vegetation thermal buffering 24 effects at lowland sites, significant drought impacts on host plant transpiration and amplified 25 leaf surface temperature, as well as reduced leaf quality linked to the seasonal advance of plant 26 phenology. 29 1 June maximum daily temperature values recorded by host plant thermal sensors at open 2 lowland microsites ranged from 22 to 42 ºC (Fig. 3a-d). Experimental studies in Pieris 3 butterflies in thermally variable environments have been conducted (i.e. treatments of short- 4 term heat stress exposure in daily cycles, mimicking natural daily variability). These treatments 5 report strong negative effects on larval growth rates and consumption rates for temperatures 6 above 39 ºC (Kingsolver, 2000; Kingsolver et al., 2006). In addition, it has been recently 7 reported that thermal conditions above 35 ºC can significantly increase egg and young larvae 8 mortality in other model species (Klockmann et al., 2016). In line with this finding, we 9 observed a significant negative effect of reduced larval size in thermal death time responses in 10 Pieris napi (F ratio = 9.67, p=0.031, Table S17), indicating a significantly increased 11 susceptibility of younger larvae to thermal stress. In the case of P. napi habitats we observed 12 that at open lowland microsites a large percentage of the maximum daily temperature records 13 (97.5%) were below 37.5 ºC and more than 90% were below 35 ºC (i.e. most of the values were 14 not surpassing the thermal threshold of 34.5 ºC estimated in TDT experiments for a daily 15 exposure of 6h to maximum temperatures). Therefore, and according to the available 16 experimental evidence in Pieris butterflies and other species, these thermal regimes should not 17 necessarily impose a strong negative impact on the survival, consumption rates and growth rates 18 of larvae if conditions of optimal host plant quality and reduced leaf drought stress were 19 simultaneously met (Kingsolver, 2000; Kingsolver et al., 2006). In line with these findings, we 20 measured reduced leaf quality and water stress conditions in host plants at open lowland 21 microsites during the summer period. We observed that in open exposed sites amplified leaf 22 surface temperatures and reduced transpiration could significantly increase host plant 23 temperature (Fig 4), surpassing the experimental thresholds estimated (i.e. TDT for TE6h and 24 TAB). Of note, previous experimental works demonstrate significant interactions in combined 25 heat stress and altered host plant quality treatments, often resulting in stronger negative impacts 26 on butterfly larval growth rates (e.g. Jones, 1982; Kingsolver, 2000). In addition, pupal mass 27 has been positively associated with fitness and total lifetime egg production in the genus Pieris 30 1 (Jones, 1982; Wiklund & Kaitala, 1995). Consequently, direct negative impacts of body size 2 reductions on population demography should not be discarded. 3 Our results also documented summer drought impacts on host plant ecophysiology. The 4 observed reductions of leaf conductances at the peak of summer drought (values < 0.2 mol/m2 s) 5 are in line with the quantitative values reported for Brassica rapa in comprehensive water stress 6 experimental treatments (Fig 5; Anissa et al., 2013, Guo et al., 2015). Under strong drought 7 stress, host plants are expected to progressively reduce leaf water content and transpiration. This 8 could potentially affect butterfly population demography because leaf-water content is known to 9 be an important factor for larval development (Soo Hoo & Fraenkel, 1966; Scriber, 1977; 10 Slansky & Feeny, 1977). Moreover, leaf transpiration and leaf water content are key characters 11 driving host plant selection by females in Pieris butterflies (Wolfson, 1980; Myers, 1985). The 12 same is true for leaf nitrogen content, which also limits larval development and is a key trait in 13 female host plant selection (Myers, 1985). Finally, leaf water and nitrogen content are generally 14 positively correlated in Brassica host plants used by Pieris species (Mattson, 1980) and are also 15 positively and significantly related to transpiration rates (Myers, 1985). On top of this, our 16 results indicated a key role of decreased June rains on long-term population declines and in 17 addition reported a significant reduction of leaf conductance in the transition from May to June 18 at lowland areas (Fig. 5). Our study also highlights the potential importance of seasonal trends 19 in leaf phenology, which in turn determine C/N content and host plant quality (Kriedeman, 20 1968; Kause et al., 1997). To our knowledge, these factors have been seldom considered as 21 contributing factors determining butterfly population vulnerability to increased drought impacts. 22 Overall this study identifies multiannual trends in summer drought as a primary driver of long- 23 term demographic declines of Pieris napi. Crucially, nearly 70% of the butterfly species in this 24 hotspot region for European butterflies are currently affected by significant population declines 25 (Stefenescu et al., 2011ab, Melero et al., 2016). Landscape changes and climatic drivers have 26 been considered as the principal candidate drivers of these widespread declines but their relative 27 role and the ecological mechanisms implied are still poorly described for most of the species. In 31 1 this context, our study clarifies the importance of summer drought as a key primary driver in P. 2 napi in the studied populations and sheds some light into some of the ecological mechanisms 3 implied (i.e. vegetation thermal buffering, phenology effects on plant quality (C/N) and changes 4 in host plant water transpiration and content). It remains to be assessed whether these 5 mechanisms could also apply to other populations of P. napi in Catalonia and to other butterfly 6 species. In this regard, it is important to bear in mind that our analyses are restricted to 7 abundant populations located in protected areas. The reported trends could possibly differ in 8 lowland and mid-land populations currently affected by increased urbanisation pressures, 9 intensified land use changes, pesticide management impacts and land abandonment (Stefanescu 10 et al 2011ab). Moreover, the results do not describe the responses of P. napi populations that 11 rely on other host plants in Catalonia (e.g. Cardamine pratensis, Arabis glabra). The host plant- 12 specific mechanisms described in the paper may non-necessarily apply to these populations. 13 Finally, our study suggests that wing and body size measures are reliable phenotypic biomarkers 14 of the geographic variability of thermal stress exposure in the studied populations, providing an 15 indirect indicator of limited habitat thermal buffering conditions for these specific populations. 16 In our field and experimental datasets, the percentage of reduction of wing size per degree 17 Celsius (as defined in Forster et al., 2012) was in the range of 1-2% [regression slope for 18 normalised experimental data: -1.56±0.25, p<0.0001; regression slope for normalised field data: 19 -1.80±0.11, p<0.0001]. This trend is consistent with the experimental slopes reported for body 20 size-temperature relationships in other temperate butterfly species exposed to similar 21 experimental thermal treatments (Forster et al. 2012; see Fig. S14 for some examples). More 22 detailed quantitative studies of the effects of thermal stress on survival and fecundity functions 23 in this species are required to estimate the critical size values associated to negative effects on 24 insect performance and the associated thermal threshold (P* and T* values, see supplementary 25 text T1 for further discussion). Moreover, we show that complementary analyses of host plant 26 dynamics are highly informative and necessary, due to the multiple ecological processes that 27 seem to be co-acting and interacting (Nygrin et al., 2008, Talloen et al., 2009). In summary, this 32 1 study indicates that phenotypic thermal biomarkers are informative as climatic stress indicators 2 but should be complemented, whenever possible, by multi-trait frameworks analysing host plant 3 ecophysiological responses and by detailed microclimatic measurements. 4 5 Acknowledgements: 6 Emili Bassols and Francesc Xavier Santaeufemia provided support with permission 7 management, scientific advice and key assistance during field work. Thanks to Ty Hedrick and 8 Heidi McLean for providing guidance and help with the matlab PhotoGUI applications 9 (Kingsolver Lab, University of North Carolina). Joel Kingsolver provided useful comments that 10 largely improved the manuscript. Thanks to Enrico Rezende and Mauro Santos for technical 11 advice in the experimental design of TDT experiments. Jordi Artola and Andreu Ubach 12 provided invaluable help collecting P. napi samples at Can Jordà and lowland sites. Melodia 13 Tamayo, Joaquim de Gispert and Andreu Ubach contributed to the experimental work. Joan 14 Llusià, Gerard Farré and Daijun Liu helped with the plant photoshynthesis measurements. 15 Thanks to Consorci per a la Protecció i Gestió dels Espais Naturals del Delta del Llobregat, Parc 16 Natural de la Zona Volcànica de la Garrotxa, and Parc Natural Aiguamolls de l’Empordà for 17 logistic support. This research was supported by VENI-NWO 863.11.021, the Spanish 18 government grants CGL2016-78093-R, CGL2013-48074-P and CGL2013-48277-P, the Catalan 19 Government project SGR 2014-274, and the European Research Council Synergy Grant ERC- 20 2013-SyG 610028 IMBALANCE-P. Additional funding was provided to CWW from the Knut 21 and Alice Wallenberg Foundation (KAW 2012.0058) and the Swedish Research Council grant 22 VR-2012-4001. 23 References: 24 Aalberg Haugen, I. M., Berger, D., & Gotthard, K. (2012). The evolution of alternative 25 developmental pathways: footprints of selection on life-history traits in a butterfly. Journal of 26 Evolutionary Biology, 25, 1377-1388. 27 33 1 Adams, M. A., & Grierson, P. F. (2001). Stable isotopes at natural abundance in terrestrial plant 2 ecology and ecophysiology: an update. Plant Biology, 3, 299–310. 3 4 Anissa, A., Chen, S., Turner, N. C., & Cowling, W. A. (2013). Genetic variation for heat 5 tolerance during the reproductive phase in Brassica rapa. Journal of Agronomy and Crop 6 Science, 199, 424-435. 7 8 Ashton, S., Gutierrez, D., & Wilson, R. J. (2009). Effects of temperature and elevation on 9 habitat use by a rare mountain butterfly: implications for species responses to climate change. 10 Ecological Entomology, 34, 437-446. 11 12 Atkinson, D. (1994). Temperature and organism size—a biological law for ectotherms. 13 Advances in Ecological Research, 25, 1–58. 14 15 Atkinson, D., & Sibly, R. M. (1997). Why are organisms usually bigger in colder 16 environments? Making sense of a life history puzzle. Trends in Ecology and Evolution, 12, 17 235–239. 18 19 Bauerfeind, S. S., & Fischer, K. (2014). Simulating climate change: temperature extremes but 20 not means diminish performance in a widespread butterfly. Population Ecology, 56, 239–250. 21 22 Bauerfeind, S. S., & Fischer, K. (2013a). Increased temperature reduces herbivore host plant 23 quality. Global Change Biology, 19, 3272–3282. 24 34 1 Bauerfeind, S. S., & Fischer, K. (2013b). Testing the plant stress hypothesis: stressed plants 2 offer better food to an insect herbivore. Entomologia Experimentalis et Applicata, 149, 148– 3 158. 4 5 Carnicer, J., Wheat, C., Vives-Ingla, M., Ubach, A., Domingo, C., Nylin, S., ... Peñuelas, J. 6 (2017). Evolutionary Responses of Invertebrates to Global Climate Change: the Role of Life- 7 History Trade-Offs and Multidecadal Climate Shifts. In S. N. Johnson & T. H. Jones (Eds.), 8 Global climate change and terrestrial invertebrates (pp. 317-348). Chichester, UK:John Wiley 9 & Sons. 10 11 Carnicer, J., Brotons, L., Stefanescu, C., & Peñuelas, J. (2012). Biogeography of species 12 richness gradients: linking adaptive traits, demography and diversification. Biological Reviews, 13 87, 457–479. 14 15 Carnicer, J., Stefanescu, C., Vila, R., Dincă, V., Font, X., & Peñuelas, J. (2013). A unified 16 framework for diversity gradients: the adaptive trait continuum. Global Ecology and 17 Biogeography, 22, 6–18. 18 19 Carnicer, J., Barbeta, A., Sperlich, D., Coll, M., & Peñuelas, J. (2013). Contrasting trait 20 syndromes in angiosperms and conifers are associated with different responses of tree growth to 21 temperature on a large scale. Frontiers in Plant Science, 4, 409. 22 23 De Frenne, P., Rodríguez-Sánchez, F., Coomes, D. A., Baeten, L., Verstraeten, G., Vellend, M., 24 ... Decocq, G. M. (2013). Microclimate moderates plant responses to macroclimate warming. 25 Proceedings of the National Academy of Sciences, 110, 18561–18565. 26 35 1 Deutsch, C. A., Tewksbury, J. J., Huey, R. B., Sheldon, K. S., Ghalambor, C. K., Haak, 2 D. C., & Martin, P. R. (2008). Impacts of climate warming on terrestrial ectotherms 3 across latitude. Proceedings of the National Academy of Sciences, 105, 6668-6672. 4 5 Domingo-Marimón, C. (2016). Contributions to the knowledge of the multitemporal spatial 6 patterns of the Iberan Peninsula droughts from a Geographic Information Service perspective. 7 PhD Thesis. Autonomous University of Barcelona. 8 9 Fischer, K., & Fiedler, K. (2000). Response of the copper butterfly Lycaena tityrus to increased 10 leaf nitrogen in natural food plants: evidence against the nitrogen limitation hypothesis. 11 Oecologia, 124, 235-241. 12 13 Forster, J., Hirst, A. G., & Atkinson, D. (2012). Warming-induced reductions in body size are 14 greater in aquatic than terrestrial species. Proceedings of the National Academy of Sciences, 15 109, 19310–19314. 16 17 García-Barros, E., Munguira, M. L., Stefanescu, C. & Vives-Moreno, A. (2013). Lepidoptera 18 Papilionoidea. Museo Nacional de Ciencias Naturales, Spain. 19 20 Grant, K., Kreyling, J., Dienstbach, L., Beierkuhnlein, C., & Jentsch, A. (2014). Water stress 21 due to increased intra-annual precipitation variability reduced forage yield but raised forage 22 quality of a temperate grassland. Agriculture Ecosystems and Environment, 186, 11–22 23 24 Guo, Y. M., Turner, N. C., Chen, S., Nelson, M. N., Siddique, K. H. M., & Cowling, W. A. 25 (2015). Genotypic variation for tolerance to transient drought during the reproductive phase of 26 Brassica rapa. Journal of Agronomy and Crop Science, 201, 267–279. 36 1 2 Hadfield, J. D. (2010). MCMC methods for multi-response generalized linear mixed models: 3 the MCMCglmm R package. Journal Statistical Software, 33, 1–22. 4 5 Jacobs, J. (2007). Ecology and Management of Whitetop (Cardaria Draba (L.) Desv.). US 6 Department of Agriculture, Natural Resources Conservation Service. 7 8 Jones, R. E., Hart, J. R., & Bull, G. D. (1982). Temperature, Size and Egg Production in the 9 Cabbage Butterfly, Pieris rapae L. Australian Journal of Zoology, 30, 223–232. 10 11 Kaitaniemi, P., Ruohomaèki, K., Ossipov, V., Haukioja, E., & Pihlaja, K. (1998). Delayed 12 induced changes in the biochemical composition of host plant leaves during an insect outbreak. 13 Oecologia, 116, 182–190. 14 15 Kause, A., Ossipov, V., Haukioja, E., Lempa, K., Hanhimäki, S., & Ossipova, S. (1999). 16 Multiplicity of biochemical factors determining quality of growing birch leaves. Oecologia, 17 120, 102–112. 18 19 Kingsolver, J. G. (2000). Feeding, growth, and the thermal environment of cabbage white 20 caterpillars, Pieris rapae L. Physiological and Biochemical Zoology, 73, 621–628. 21 22 Kingsolver, J. G., Shlichta, J. G., Ragland, G. J., Massie, K. R. (2006). Thermal reaction norms 23 for caterpillar growth depend on diet. Evolutionary Ecology Research, 8, 703–715. 24 25 Kingsolver, J. G. (2009). The well-temperatured biologist. American Naturalist, 174, 755–768. 26 37 1 Kirk, J. T. O. (1994). Light and photosynthesis in aquatic ecosystems. Cambridge, UK: 2 Cambridge University Press. 3 4 Klockmann, M., Günter, F., & Fischer, K. (2016) Heat resistance throughout ontogeny: body 5 size constrains thermal tolerance. Global Change Biology, early view. 6 7 Kriedemann, P. E. (1968). Photosynthesis in vine leaves as a function of light intensity, 8 temperature, and leaf age. Vitis, 7, 213–220. 9 10 Marshall, J. D., Brooks, J. R., Lajtha, K. (2007). Sources of variation in the stable isotopic 11 composition of plants. Stable Isotopes in Ecology and Environmental Science, 2, 22–60. 12 13 Mattson, W. J. (1980). Herbivory in relation to plant nitrogen content. Annual Review of 14 Ecology, Evolution and Systematics, 11, 119–161. 15 16 Melero, Y., Stefanescu, C., & Pino, J. (2016). General declines in Mediterranean butterflies 17 over the last two decades are modulated by species traits. Biological Conservation, 201, 336– 18 342. 19 20 Merrill, R. M., Gutiérrez, D., Lewis, O.T., Gutiérrez, J., Díez, S. B., Wilson, R. J. (2008). 21 Combined effects of climate and biotic interactions on the elevational range of a phytophagous 22 insect. Journal of Animal Ecology, 77, 145–155. 23 24 Mulligan, G. A., & Findlay, J. N. (1974). The biology of Canadian weeds. Canadian Journal of 25 Plant Science, 54, 149–160. 26 38 1 Munns, R., James, R. A., Sirault, X. R., Furbank, R. T., & Jones, H. G. (2010). New 2 phenotyping methods for screening wheat and barley for beneficial responses to water deficit. 3 Journal of Experimental Botany, 61, 3499–3507. 4 5 Myers, J. H. (1985). Effect of physiological condition of the host plant on the ovipositional 6 choice of the cabbage white butterfly, Pieris rapae. Journal of Animal Ecology, 54, 193–204. 7 8 Nail, K. R., Batalden, R. V., & Oberhauser, K. S. (2015). What's too hot and what's too cold? 9 In: K. S. Oberhauser, K. R. Nail, & S. Altizer (Eds.), Monarchs in a changing world: biology 10 and conservation of an iconic butterfly (pp. 99–108). New York, USA: Cornell University 11 Press. 12 13 Nieto-Sánchez, S., Gutiérrez, D., & Wilson R. J. (2015). Long-term change and spatial variation 14 in butterfly communities over an elevational gradient: driven by climate, buffered by habitat. 15 Diversity and Distributions, 21, 950–961. 16 17 Nygren, G. H., Bergström, A., & Nylin, S. (2008). Latitudinal body size clines in the butterfly 18 Polyommatus icarus are shaped by gene-environment interactions. Journal of Insect Science, 8, 19 1-13. 20 21 Oliver, T. H., Stefanescu, C., Páramo, F., Brereton, T., & Roy, D. B. (2014). Latitudinal 22 gradients in butterfly population variability are influenced by landscape heterogeneity. 23 Ecography, 37, 863-871. 24 39 1 Oliver, T. H., Marshall, H. H., Morecroft, M. D., Brereton, T., Prudhomme, C., & Huntingford, 2 C. (2015). Interacting effects of climate change and habitat fragmentation on drought-sensitive 3 butterflies. Nature Climate Change, 5, 941-945. 4 5 Parmesan, C., Root, T. L., & Willig, M. R. (2000). Impacts of extreme weather and climate on 6 terrestrial biota. Bulletin of the American Meteorological Society, 81, 443–450. 7 8 Pateman, R. M., Thomas, C. D., Hayward, S. A., & Hill, J. K. (2016). Macro and microclimatic 9 interactions can drive variation in species' habitat associations. Global Change Biology, 22, 10 556–566. 11 12 Peñuelas J, & Filella, I. (2001). Responses to a warming world. Science, 294, 793–795. 13 14 Pollard, E., Yates, T. J. (1994). Monitoring butterflies for ecology and conservation: the British 15 butterfly monitoring scheme. London, UK: Springer Science & Business Media. 16 17 Pons, X. (2002). MiraMon. Geographic Information System and Remote Sensing software. 18 CREAF. Bellaterra, Spain. ISBN: 84-931323-5-7. 19 20 Rezende, E. L., Castañeda, L. E., & Santos, M. (2014). Tolerance landscapes in thermal 21 ecology. Functional Ecology, 28, 799-809. 22 23 SAS Institute Inc. (2012). JMP 10. Cary, NC, SAS Institute Inc. 24 40 1 Scriber, J. M. (1977). Limiting effects of low leaf-water content on the nitrogen utilization, 2 energy budget and larval growth of Hyalophora cecropia (Lepidoptera: Saturnidae). Oecologia, 3 28, 264–287. 4 5 Scriber, J. M., Slansky, F. (1981). The nutritional ecology of immature insects. Annual Review 6 of Entomology, 26, 183–211. 7 8 Schmucki, R., Pe'Er, G., Roy, D. B., Stefanescu, C., Van Swaay, C. A., Oliver, T. H., ... 9 Musche, M. (2016). A regionally informed abundance index for supporting integrative analyses 10 across butterfly monitoring schemes. Journal of Applied Ecology, 53, 501–510. 11 12 Sheridan, J. A., & Bickford, D. (2011). Shrinking body size as an ecological response to climate 13 change. Nature Climate Change, 1, 401–406. 14 15 Slansky, F. J., & Feeny, P. (1977). Stabilization of the rate of nitrogen accumulation by larvae 16 of the cabbage butterfly on wild and cultivated food plants. Ecological Monographs, 47, 209– 17 228. 18 19 Soo Hoo, C. F., & Fraenkel, G. (1966). The consumption, digestion and utilization of food 20 plants by a polyphagous insect Prodenia eridania (Cramer). Journal of Insect Physiology, 12, 21 711–730. 22 23 Stefanescu, C., Herrando S., & Páramo F. (2004). Butterfly species richness in the north-west 24 Mediterranean Basin: The role of natural and human-induced factors. Journal of Biogeography, 25 31, 905–915. 26 41 1 Stefanescu, C., Torre, I., Jubany, J., & Páramo, F. (2011a). Recent trends in butterfly 2 populations from north-east Spain and Andorra in the light of habitat and climate change. 3 Journal of Insect Conservation, 15, 83–93. 4 5 Stefanescu, C., Carnicer, J., & Penuelas, J. (2011b). Determinants of species richness in 6 generalist and specialist Mediterranean butterflies: the negative synergistic forces of climate and 7 habitat change. Ecography, 34, 353–363. 8 9 Suggitt, A. J., Stefanescu, C., Páramo, F., Oliver, T., Anderson, B. J., Hill, J. K., ... Thomas, C. 10 D. (2012). Habitat associations of species show consistent but weak responses to climate. 11 Biology Letters, 8, 590–593. 12 13 Suggitt, A. J., Wilson, R. J., August, T. A., Fox, R., Isaac, N. J., Macgregor, N. A., ... Maclean, 14 I. M. (2015). Microclimate affects landscape level persistence in the British Lepidoptera. 15 Journal of Insect Conservation, 19, 237–253. 16 17 Sunday, J. M., Bates, A. E., Kearney, M. R., Colwell, R. K., Dulvy, N. K., Longino, J. T., & 18 Huey, R. B. (2014). Thermal-safety margins and the necessity of thermoregulatory behavior 19 across latitude and elevation. Proceedings of the National Academy of Sciences, 111, 5610– 20 5615. 21 22 Talloen, W., Van Dongen, S., Van Dyck, H., & Lens, L. (2009). Environmental stress and 23 quantitative genetic variation in butterfly wing characteristics. Evolutionary ecology, 23, 473- 24 485. 25 42 1 Valim, J. O. S., Teixeira, N. C., Santos, N. A., Oliveira, M. G. A., & Campos, W. G. (2016). 2 Drought-induced acclimatization of a fast-growing plant decreases insect performance in leaf- 3 chewing and sap-sucking guilds. Arthropod-Plant Interactions, 10, 351–363. 4 5 Wahlberg, N., Rota, J., Braby, M. F., Pierce, N. E., & Wheat, C. W. (2014). Revised 6 systematics and higher classification of pierid butterflies (Lepidoptera: Pieridae) based on 7 molecular data. Zoologica Scripta, 43, 641-650. 8 9 Wiklund, C., & Kaitala, A. (1995). Sexual selection for large male size in a polyandrous 10 butterfly: the effect of body size on male versus female reproductive success in Pieris napi. 11 Behavioral Ecology, 6, 6–13. 12 13 Wilson, R. J., Gutiérrez, D., Gutiérrez, J., Martínez, D., Agudo, R., & Monserrat, V. J. (2005). 14 Changes to the elevational limits and extent of species ranges associated with climate change. 15 Ecology Letters, 8, 1138–1146. 16 17 Wilson, R. J., Gutierrez, D., Gutierrez, J., & Monserrat, V. J. (2007). An elevational shift in 18 butterfly species richness and composition accompanying recent climate change. Global Change 19 Biology, 13, 1873–1887. 20 21 Wolfson, J. L. (1980). Oviposition response of Pieris rapae to environmentally induced 22 variation in Brassica nigra. Entomologia Experimentalis et Applicata, 27, 223–232. 23 24 Zografou, K., Kati, V., Grill, A., Wilson, R. J., Tzirkalli, E., Pamperis, L. N., & Halley, J. M. 25 (2014). Signals of climate change in butterfly communities in a Mediterranean protected area. 26 PloS one, 9, e87245. 43 1 44