Received: 1 July 2022 Revised: 4 November 2022 Accepted: 14 November 2022 DOI: 10.1002/ecm.1561 ARTICLE Interspecific differences in microhabitat use expose insects to contrasting thermal mortality Maria Vives-Ingla 1,2 | Javier Sala-Garcia 1,3 | Constantí Stefanescu 1,4 | Armand Casado -Tortosa 5 | Meritxell Garcia 5 | Josep Peñuelas 1,6 | Jofre Carnicer 2,5,7 1 CREAF, Barcelona, Spain 2 Abstract Universitat Autònoma de Barcelona, Barcelona, Spain Ecotones linking open and forested habitats contain multiple microhabitats 3 Departamento de Biodiversidad y with varying vegetal structures and microclimatic regimes. Ecotones host Biología Evolutiva, Museo Nacional de many insect species whose development is intimately linked to the microcli- Ciencias Naturales, CSIC, Madrid, Spain 4 matic conditions where they grow (e.g., the leaves of their host plants and the Museu de Ciències Naturals de Granollers, Granollers, Spain surrounding air). Yet microclimatic heterogeneity at these fine scales and its 5 Departament de Biologia Evolutiva, effects on insects remain poorly quantified for most species. Here we studied Ecologia i Ciències Ambientals, Facultat how interspecific differences in the use of microhabitats across ecotones lead de Biologia, Universitat de Barcelona, Barcelona, Spain to contrasting thermal exposure and survival costs between two closely-related 6 CSIC, Global Ecology butterflies (Pieris napi and P. rapae). We first assessed whether butterflies CREAF-CSIC-UAB, Cerdanyola del selected different microhabitats to oviposit and quantified the thermal condi- Vallès, Spain tions at the microhabitat and foliar scales. We also assessed concurrent 7 Institut de Recerca de la Biodiversitat changes in the quality and availability of host plants. Finally, we quantified (IRBio), Universitat de Barcelona, Barcelona, Spain larval time of death under different experimental temperatures (thermal death time [TDT] curves) to predict their thermal mortality considering both the Correspondence Jofre Carnicer intensity and the duration of the microclimatic heat challenges in the field. Email: jofre.carnicer@ub.edu We identified six processes determining larval thermal exposure at fine scales associated with butterfly oviposition behavior, canopy shading, and heat and Funding information Agència de Gestio d’Ajuts Universitaris i water fluxes at the soil and foliar levels. Leaves in open microhabitats could de Recerca, Grant/Award Number: reach temperatures 3–10 C warmer than the surrounding air while more SGR-2017-1005; Ministerio de Asuntos closed microhabitats presented more buffered and homogeneous temperatures. Econ omicos y Transformaci on Digital, Gobierno de España, Grant/Award Interspecific differences in microhabitat use matched the TDT curves and the Numbers: CGL2013-48074-P, thermal mortality in the field. Open microhabitats posed acute heat challenges CGL2016-78093-R; Ministerio de Ciencia e that were better withstood by the thermotolerant butterfly, P. rapae, where the Innovacion, Grant/Award Numbers: FPU17/05869, PID2020-117636GB-C21; species mainly laid their eggs. Despite being more thermosensitive, P. napi Nederlandse Organisatie voor was predicted to present higher survivals than P. rapae due to the thermal Wetenschappelijk Onderzoek, Grant/Award Number: 863.11.021 buffering provided by their selected microhabitats. However, its offspring could be more vulnerable to host-plant scarcity during summer drought Handling Editor: Matthew L. Forister periods. Overall, the different interaction of the butterflies with microclimatic This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2022 The Authors. Ecological Monographs published by Wiley Periodicals LLC on behalf of The Ecological Society of America. Ecological Monographs. 2023;93:e1561. https://onlinelibrary.wiley.com/r/ecm 1 of 19 https://doi.org/10.1002/ecm.1561 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 2 of 19 VIVES-INGLA ET AL. and host-plant variation emerging at fine scales and their different thermal sensitivity posed them contrasting heat and resource challenges. Our results contribute to setting a new framework that predicts insect vulnerability to cli- mate change based on their thermal sensitivity and the intensity, duration, and accumulation of their heat exposure. KEYWORDS microclimates, microhabitat, Pieris, plant–insect interactions, thermal adaptations, thermal mortality, thermal tolerance INTRODUCTION i.e., microhabitat), while the microclimate measured at the plant surfaces might be more important for eggs and small Vegetation cover locally modifies climatic conditions larvae (Kingsolver et al., 2011; Pincebourde et al., 2021; and generates a microclimatic regime that deviates from Potter et al., 2009; Woods, 2010, 2013). Microclimatic het- open, free-air, and standardized measurements (i.e., erogeneity at these fine scales can be particularly high and macroclimate; Geiger, 1950; Stoutjesdijk & Barkman, 2014). comparable to that recorded at macroclimatic scales, but The absorption and reflection of solar radiation and the remains poorly quantified for most species (Pincebourde evapotranspirative cooling of forest canopies buffer et al., 2016; Pincebourde & Woods, 2020). macroclimatic temperatures and reduce thermal variation The selection of microhabitat and host plant during in the understory (Bramer et al., 2018; De Frenne butterfly oviposition will affect the growing environment et al., 2021; Zellweger et al., 2020). In contrast, in open of the offspring (Figure 1a; Courtney, 1986; Doak areas with short and sparse vegetation, temperatures near et al., 2006; Forsberg, 1987; Gibbs & Van Dyck, 2009). On the ground can be more extreme than those recorded at the one hand, we expect that shadier microhabitats will 2-m and shady conditions (i.e., thermal amplification; offer buffered microclimates, with dampened thermal Carnicer et al., 2021; Stoutjesdijk & Barkman, 2014; Woods variability and extreme values, while they are amplified et al., 2015). Narrow ecotones that link open and forested in open areas exposed to direct radiation (Figure 1a habitats (hereafter termed open–closed ecotones) generate A–C). On the other hand, we also expect that microhabi- multiple microhabitats with varying vegetal structures, tat conditions could induce plastic shifts on many traits potentially exposing the organisms they harbor to that define host plant quality, depending on their shade contrasting microclimatic conditions. Here we studied how tolerance (Figure 1a D–E; Poorter et al., 2019; Scriber & interspecific differences in microhabitat use across ecotones Slansky, 1981). Variations in host plant quality can inter- shape the thermal exposure and the associated costs to sur- act with the microclimate experienced by the feeding lar- vival of two closely-related butterflies. vae and influence their development in complex ways Microclimatic variation determines the thermal experi- (Clissold & Simpson, 2015). ence of the organisms and has important effects on Butterfly responses to microclimatic exposure will their thermal adaptations and performance (Franken depend on both their thermal exposure and thermal sen- et al., 2018; Kaiser et al., 2016; Kaspari et al., 2015; sitivity (Carnicer et al., 2017). We hypothesize that the Pincebourde & Casas, 2019; Woods et al., 2022). The scale two studied butterflies will show diverging thermal adap- at which microclimatic measurements are relevant for tive strategies according to their microhabitat preferences organisms depends on their body size and mobility and the associated thermal regime (Figure 1a F–G). The (Kingsolver et al., 2011; Pincebourde et al., 2021; probability of surviving heat stress depends on both its Pincebourde & Woods, 2020). For butterflies, microclimatic intensity and its duration (Rezende et al., 2014). This variability at the landscape scale is pertinent in the adult relationship can be synthesized by thermal death time stage, as they can easily fly and sample between alternative (TDT) curves, where the critical temperature (Tko) that habitats (e.g., woodland vs. grassland; Suggitt et al., 2011, an organism can tolerate linearly decreases with the loga- 2012). In contrast, the area that eggs and larvae can use is rithm of exposure time (t) following this equation: much smaller (Courtney, 1986); thus, these less mobile stages will be more likely affected by the microclimatic var- T ko ¼ CT max  z log 10 t: ð1Þ iability inside their habitats. Larval stages may be more likely influenced by the microclimatic conditions of the air Therefore, thermal tolerance to a heat challenge not only bath surrounding them (at less than ~1 m of distance, depends on the upper critical thermal limit (CTmax) but 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ECOLOGICAL MONOGRAPHS 3 of 19 (a) Open–closed ecotone Microhabitat and A host-plant selection Thermosensitive Thermotolerant B Thermal F C buffering G Thermal amplification Shade-tolerant Light-demanding host plant host plant D E Thermal exposure and microclimatic regimes Gradients Host-plant condition Larval thermal tolerance and thermal mortality (b) y-intercept = tmax= CTmax /z CT max 1 t > tmax’ log 10 t = T z z Survival time (log10 t) max Thermosensitive strategy slope = 1/z z < z’ Thermotolerant strategy (’) x-intercept = CTmax T* CTmax < CTmax’ Thermal Chronic Acute Temperature (T) regime stress stress Moderate and Extreme and less variable more variable temperatures temperatures F I G U R E 1 (a) Narrow ecotones that link open and forested habitats generate multiple microhabitats with distinct microclimatic regimes. Insects with contrasting thermal tolerances may select different microhabitats to oviposit, leading to large interspecific differences in the thermal exposure of eggs and larvae. Moreover, host plants with varying shade tolerances may show different distributions, traits, and conditions across ecotones, synergistically affecting larval performance. Overall, concurrent gradients of both biotic and abiotic processes are produced across ecotones. (b) Thermotolerant and thermosensitive strategies are expected to evolve under different thermal regimes. The two strategies are characterized by different thermal death time (TDT) curves, reflecting an evolutionary trade-off between survival capacity at acute, extreme stresses (thermotolerant) and chronic, less intense conditions (thermosensitive). Both species exhibit equal survival times at T*, which represents the thermal threshold between these two alternative thermal strategies. Silhouettes used in this figure were obtained from rawpixel.com and phylopic.org. 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 4 of 19 VIVES-INGLA ET AL. also on the thermal sensitivity (z), which describes the long-term monitoring network that provides data on required increase of temperature to decrease time to butterfly abundance at a weekly resolution since 1994 death by one order of magnitude (Figure 1b). The analy- (i.e., CBMS; Pollard & Yates, 1993; Stefanescu, 2000). They sis of the TDT curves of 56 species of insects, bivalves, contain abundant populations of the green-veined white and fishes has pointed out that CTmax and z are positively (P. napi) and the small white (P. rapae) that have been and tightly associated (Rezende et al., 2014). Organisms intensively studied since 2012 (Carnicer et al., 2019; with a high CTmax and z are more capable of surviving Vives-Ingla et al., 2020). Site 1 is in a mid-elevation area extreme temperatures, but their critical temperature (Tko) (539 m above sea level (asl); Can Jordà, La Garrotxa rapidly decreases with exposure time (Figure 1b, yellow). Volcanic Zone Natural Park) with a heterogeneous land- In contrast, low z values allow longer survival times at scape of evergreen and deciduous woodlands, meadows, less intense but still stressful temperatures to the detri- and natural ponds. Site 2 is in a coastal wetland (2 m asl; ment of CTmax (Figure 1b, blue). As a consequence of this El Cortalet, Aiguamolls de l’Empordà Natural Park) trade-off, we expect that larvae growing in open micro- surrounded by riparian deciduous forests, bush and bram- habitats may show a thermotolerant strategy (high CTmax ble thickets, reed beds, and irrigated cropland. The land- and z) to cope with the acute stresses of their scape mosaic of both sites generates spatial gradients of highly-variable thermal regime (extreme temperatures vegetation cover between open and closed habitats and for short times). In the same line, we would expect the their transition zones (open–closed ecotones). opposite, thermosensitive strategy for species selecting P. napi and P. rapae are partially syntopic species shadier microhabitats with less intense temperatures but (i.e., they share some of the microhabitats for breeding). longer exposures (i.e., chronic stress). They co-occur in a wide range of habitats from the sea To assess whether microhabitat variability in level to alpine regions, but P. napi is usually linked to open–closed ecotones generates interspecific differences shaded, humid sites while P. rapae is more common in in their thermal exposure and the associated costs to dry, open areas. Both butterflies lay individual eggs on survival, we studied two model species of butterfly: Brassicaceae species, such as A. petiolata and L. draba. Pieris napi L. 1758 and P. rapae L. 1758. These species However, P. rapae is a more generalist species and uses a were studied in two Mediterranean sites, which harbor greater diversity of host plants (Ohsaki, 1979). populations of two host plants with contrasting toler- The dominant host plant in the mid-elevation site is ances to shade, Alliaria petiolata (shade-tolerant) and Alliaria petiolata. It is a biennial herb common in damp, Lepidium draba (light-demanding). We first assessed shaded soils at the edges of deciduous and riverine for- whether the two butterfly species selected different ests. It can grow in highly contrasted environmental con- microhabitats from the open–closed ecotones to oviposit ditions, exhibiting considerable plasticity in different and quantified the thermal exposure of their offspring habitats (Cavers et al., 1979). Seedlings emerge during both at the microhabitat and foliar scales. We also spring and early summer and persist as rosettes through- assessed concurrent changes in the nutritional quality out the first year, until the next growing season when and condition of host plants across ecotones. Then, we inflorescences are initiated. Lepidium draba is the domi- conducted ecophysiological assays of heat tolerance nant host plant in the lowland site. It is a perennial, rhi- with the larvae of both species to estimate their TDT zomatous herb and can be found in open areas and field curves and determine their thermal sensitivity. Based margins (de Bolo s & Vigo, 1990). Its extensive, on the experimental quantification of larval survival at multibranched rhizomes are notably capable of produc- different temperatures, we finally applied a dynamic ing many new shoots, which can develop into large model to predict thermal mortality in field microclimatic monocultural stands (Francis & Warwick, 2008). conditions. Oviposition behavior MATERIALS AND METHODS We assessed whether the differences in broad habitat Study system preferences between the two butterflies (P. napi for humid areas and P. rapae for open habitats) led to differ- The thermal exposure of the two species, its impacts, and ent microhabitat selection for oviposition across the concurrent variation in their host plants were studied open–closed ecotones. We tested this hypothesis by carry- in two protected areas of the north-eastern Iberian ing out censuses of behavior at the two study sites. Peninsula, 50 km from each other (Appendix S1: Females were followed and their behaviors were recorded Figures S1 and S2). The sites of study are part of a for periods of 45 min. The censuses fully covered the 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ECOLOGICAL MONOGRAPHS 5 of 19 entire daily period of flight activity, between 9:00 and Every monitoring day, we selected at least 16 host 19:00, and were conducted during summer. Oviposition plants for each cohort (4 individuals  4 microhabitat was considered to occur when females that landed on a types) ensuring that plants were randomly chosen, with- leaf were observed to curl their abdomen and remain in out repetition to avoid pseudoreplication. Microclimatic this position for at least 3 s. Egg-laying was visually con- variables included several temperatures at the soil, foliar firmed in most of the cases. Species, hour, and microhab- and air level, and soil moisture. We used a penetration itat type (open, closed, or intermediate) were noted. thermometer (HI98509, Hanna Instruments Ltd, Spain) Censuses were simultaneously performed in the various to measure soil temperature at a depth of 10 cm. Soil sur- microhabitat types, balancing the time spent in each type. face temperature, air temperature above the host plant, The temperature of the leaves where eggs were laid and and foliar surface temperature were measured using a their position (upper vs. underside of leaves) were also thermocouple (see device and measurement details in the recorded when possible using a wire K-type thermocou- section Oviposition behavior). A minimum of three repli- ple probe (SC-TT-KI-30-1M, Omega Engineering Ltd, cates were taken for each measurement. The tempera- UK) attached to a hand-held thermocouple thermometer tures were measured between 10:00 and 16:00, and the (HH503, Omega Engineering Ltd, UK, and HI935005N, time, wind, and radiation conditions were recorded. We Hanna Instruments Ltd, Spain) immediately after the measured soil surface temperature near the host plants, female left the plant. replicating it in spots exposed to direct solar radiation and in the shade; air temperature, immediately above the host plant at a height of 1 m; and foliar temperature, on Microclimatic and host-plant variation the upper and underside of the leaves. The volumetric water content of the soil (% by volume) was measured at We assessed microclimatic conditions at both microhabitat three points near each plant using a DELTA-T SM150 and host-plant scales. For the measurements at the micro- (Delta-T Devices Ltd, UK) soil-moisture sensor. habitat scale, we installed 15 standalone data loggers The traits measured in host plants included phenology, (EL-USB-2-LCD, Lascar Electronics, UK) in different stem length, and foliar dimensions, density, and chloro- microsites harboring host plants and where ovipositing phyll and water contents. These traits were selected females had been detected. Seven of the sensors were because they have been associated with oviposition behav- installed in the mid-elevation site, and the other eight in ior and offspring performance (Awmack & Leather, 2002; the lowland site (Appendix S1: Tables S1 and S2). The sen- Gibbs & Van Dyck, 2009; Stefanescu et al., 2006; sors were programmed to record temperature ( C) and rel- Wolfson, 1980) and are simple to measure. The foliar chlo- ative humidity (%) at hourly resolution and were placed rophyll content is also an indicator of plant nutritional 25 cm above the soil protected from direct solar radiation. condition, photosynthetic capacity, and developmental To quantify thermal exposure and trait variability of stage (Curran et al., 1990; Everiit et al., 1985). We assessed host plants, we monitored cohorts of 242 individuals of plant phenological status by classifying the individuals into A. petiolata and 362 individuals of L. draba distributed in one of four phenological stages: spring early vegetative, four representative categories of microhabitats found reproductive, senescent, and summer late vegetative. A across the ecotones (open, semi-open, semi-closed, and representative basal, medial, and apical leaf was chosen for closed). Each microhabitat category was assigned based on each plant, and its state (green or senescent) was recorded. detailed measurements of the dynamics of the canopy and Chlorophyll content was estimated as the mean of three the ground cover by herbaceous plants (see Appendix S1: measurements from a SPAD-502 chlorophyll meter Section S1 and Figures S2 and S3 for further details). (Konica Minolta Sensing, Spain). Finally, leaves were sev- Canopy closure was measured by visual inspection in the ered and immediately weighed (fresh weight, FW) using a vertical and the four cardinal directions. The herbaceous field digital scale (PJS020, PESOLA Präzisionswaagen AG, layer was characterized both as herbaceous ground cover Switzerland), and then oven-dried in the laboratory at and as mean herb height using the point-intercept 60 C for 2 days to a stable weight (dry weight, DW). Foliar method. Closed and semi-closed microhabitats presented a water content was defined as (FW  DW)/DW. Foliar mean canopy closure higher than 50%, while we defined density was calculated as the ratio between DW and semi-open and open microhabitats by their herbaceous foliar length. When host plants were mature, we also layer, as trees and shrubs were less common there. We counted the number of fruits per plant (siliques for monitored the host plants in 2017, from March to A. petiolata and silicules for L. draba), as a proxy of plant September, and measured the same 14 microclimatic, phe- reproductive performance between microhabitats and nological, morphological, and physiological traits every shade tolerance. A minimum of seven individuals were 15 days (Appendix S1: Table S2). sampled for each microhabitat type. 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 6 of 19 VIVES-INGLA ET AL. Ecophysiological assays of heat tolerance were added as fixed factors, and site, date, and period of the day of the census were treated as categorical, random We implemented a static heat tolerance experiment using factors. The logarithm of census duration was added as an larvae of P. napi and P. rapae to determine whether they offset term (i.e., a linear predictor without an estimated differ in their thermal strategies. If P. napi oviposits in regression parameter to account for differences in sam- more closed and buffered microhabitats than does pling effort between censuses; Zuur et al., 2009). We P. rapae, as we hypothesized, the larvae of the former repeated the same modeling procedure using the number would exhibit a thermosensitive strategy, with lower of ovipositing females per census as the response variable values of z and CTmax. Females from both locations and (model 2). species were captured and their offspring were reared in The seasonal dynamics of microclimatic conditions growing chambers at 22 C 13L:11D, with fresh and abun- and host-plant traits were assessed by regressing LOESS dant host plants (L. draba and A. petiolata). The experi- models against ordinal day. To examine the spatial vari- ment was conducted on 210 larvae from 20 family lines ation of these variables across the open–closed ecotones, (Appendix S1: Table S3). Before the application of the an analysis of variance (ANOVA) testing for microhabi- thermal treatment, larvae were acclimated for 1 h at con- tat type was applied followed by a post hoc Tukey hon- stant 22 C and deprived of food. We recorded the larval est significant difference (HSD) test calculated using the initial weight (g) and subsequently placed the larvae in emmeans package (Lenth, 2020). The analyses were individual plastic vials (diameter: 3.5 cm, height: 7 cm) performed for the entire sampling period and specific that were submerged in a water bath programmed at a phenological stages and seasons. The characterization of constant temperature (i.e., 40, 42, or 44 C, depending on the thermal regimes in different microhabitats and the treatment). These temperatures can be recorded in scales considered the absolute records at different levels the field and are known to be stressful for both species (air, soil, leaves), thermal tendencies relative to (Kingsolver, 2000; von Schmalensee et al., 2021). The sta- macroclimatic conditions, and measures of dispersion tus of the larvae (alive or dead) was checked at regular (i.e., standard deviation and skewness). Daily records time intervals (once every 30 min for the assays at 40 C; of macroclimatic temperatures were obtained from every 20 min at 42 C; and every 10 min at 44 C), trading two meteorological stations near the study sites off between accurate detection of time to death and (Appendix S1: Figure S1). Microhabitat thermal offset potential thermal fluctuations associated with larval sta- was calculated as the difference in daily mean tempera- tus checking. The air temperature inside the plastic vial tures between the data logger and the standardized was continuously recorded using a data logger with a weather station measurements. We defined foliar ther- 20-s resolution to have a more accurate estimate of the mal offset as the difference between the upper side foliar thermal exposure of the larvae and its temporal fluctua- temperature and the synchronic air temperature above tions during the treatment. The average temperature the host plant at 1 m height (also termed “thermal recorded with the data logger was used in the subsequent excess” in some studies; see De Frenne et al., 2021; analyses, rather than the fixed, programmed temperature Pincebourde & Woods, 2012). The daily standard devia- in the water bath (which was considered a less accurate tion of temperatures recorded with the data logger was proxy). used as an indicator of the microclimatic temporal vari- ability at the microhabitat scale, while we used the daily SD of foliar temperatures in the same microhabitat as a Statistical analysis measure of spatial thermal heterogeneity at finer scales. To assess whether the two butterflies differ in their All data were analyzed using R 3.6.1 (R Core Team, 2019). thermal strategy, we fitted a linear model of the logarithm We applied generalized linear mixed models (GLMMs) to of the time to larval death against the mean temperature determine whether the two butterfly species selected dif- recorded with the data logger during the static treatments ferent microhabitats to oviposit. The final data set included for each species (i.e., TDT curve). We estimated z and a total of 7217 min of census, 139 ovipositions, and 43 CTmax from the regressed equations and we then tested ovipositing females (Appendix S1: Table S4). The number whether the two species presented different slopes by of ovipositions observed per female was used as the fitting an analysis of covariance (ANCOVA) of the effect of response variable (model 1). We fitted the model using the temperature, species, and their interaction on log10 of larval glmer function of the lme4 package (Bates et al., 2015) by survival time. Additional general mixed linear models were maximum likelihood (Laplace approximation) and a also fitted for each species on log10 of time to larval death Poisson error distribution with a log link function. The with thermal treatment ( C), larval weight (g), and site as type of microhabitat, the species, and their interaction fixed factors and family within site as a random factor. 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ECOLOGICAL MONOGRAPHS 7 of 19 Predictions of thermal mortality in the (i.e., mostly in summer) to summarize the peaks of thermal field mortality for each species in each microhabitat. Finally, we assessed how thermal mortality would change if we con- The semilogarithmic link between knockdown times and sider that larvae can avoid acute thermal stress on the temperature during a heat challenge (i.e., TDT curve) can leaves by actively moving to shadier parts of the plant be estimated experimentally for different levels of mortal- (thermal avoidance behavior, TAB; Carnicer et al., 2019). ity (e.g., the time where 100% or 50% of larvae were dead), This was done by truncating daily thermal profiles at fixed resulting in parallel lines with the same slope (z). From TAB thresholds (none, 35, 37.5, 40, 42.5, and 45 C) and these parallel curves, we can define how mortality rate repeating the whole numeric procedure with the truncated (and, hence, survival probability) changes with exposure thermal profile. These thresholds were obtained from pre- time for any constant temperature (Rezende et al., 2020b). vious experimental observations of thermal avoidance We obtain in this way temperature-specific survival behaviors for the genus Pieris (Carnicer et al., 2019). We curves, which are mathematically related and collectively performed 100 bootstrap replicates per threshold to have define the thermal tolerance landscape of the species an estimate of the uncertainty associated with the (Rezende et al., 2014, 2020a, 2020b). In the field, organ- predictions. isms are exposed to variable temperatures and their sur- vival responses are the result of shifts between temperature-specific curves. Therefore, if we know how RESULTS temperature changes during a specific period, we can pre- dict the survival probability under field conditions by sum- Interspecific differences in microhabitat ming up the infinitesimal changes in the survival rate that selection occurred during that period (Rezende et al., 2020b). Here we adapted the methods and scripts developed by Results of the GLMM analyses applied on the number of Rezende et al. (2020a, 2020b) to numerically predict the ovipositions per female (model 1) and the number of daily thermal mortality throughout the year for the larvae ovipositing females (model 2) were very similar (Table 1). of both Pieris species in the different microhabitats. Both P. napi and P. rapae distributed their eggs unequally Microclimatic field conditions during the whole year were across the open–closed ecotones (model 1: microhabitat obtained from the thermal records of the data loggers. For χ2 = 13.1, df = 2, p = 0.0015). However, the microhabitats each sensor, we calculated the mean thermal series of all selected by females differed between the two butterflies the recorded years (Appendix S1: Table S1) and estimated (model 2: microhabitat  species χ2 = 21.2, df = 2, the thermal profile at 1-min resolution of each day by p < 0.0001). P. napi preferentially selected host plants from non-linear interpolation between consecutive hours the microhabitats with intermediate vegetation covering (Rezende et al., 2020b). Daily thermal mortality was (which include the semi-open and semi-closed microhabi- predicted from March to September (both included), to tats, Figure 2a). Concretely, predictions from model 1 indi- capture the period of higher thermal stress (i.e., summer). cated that the number of eggs laid in these microhabitats We predicted larval survival in the dynamical field condi- was six and three times higher than in closed and open tions based on the TDT curves, by bootstrapping the time microhabitats, respectively. In sharp contrast, P. rapae of death of 35 larvae per species in each of the three exper- mainly laid eggs in the open microhabitats (in model 1, open imental treatments (i.e., 40, 42, or 44 C, total sample microhabitats received 99% of P. rapae ovipositions). Closed size = 210). Bootstrapped data were then used to estimate microhabitats were rarely selected by either species, but a the tolerance landscape of each species with a resolution few ovipositions of P. napi were observed. The oviposition of 0.001 of survival probability (see Rezende et al., 2020b pattern did not present relevant differences between sites, for further details). From this relationship between sur- as the variance estimated for this random factor was much vival and time at constant temperatures, we calculated the lower than the other effects in both models (model 1: daily mortality curve at a 1-min resolution for each spe- sdsite = 0.02, sddate = 1.47, sdperiod = 1.77; model 2: cies, day, and microclimatic sensor by numerical approxi- sdsite = 0.23, sddate = 0.7, sdperiod = 1.26). mation (Rezende et al., 2020b). We defined daily thermal mortality as the maximum mortality (1—survival probability) of the day and calcu- Interspecific differences in thermal lated, with these daily values, the cumulative mortality exposure during development (using a rolling window of 30 days from March to September). We then averaged cumulative Foliar temperatures recorded during oviposition differed mortality for the periods where it exceeded 0.01 between both species (with thermal differences in the 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 8 of 19 VIVES-INGLA ET AL. TABLE 1 Generalized linear mixed models of the number of ovipositions per female and the number of ovipositing females per census. Model Response variable Fixed effect χ2 df p 1 (R 2 GLMM = 0.23) Ovipositions per female Microhabitat type 13.1 2 0.0015 Species 0.3 1 0.5576 Microhabitat type  species 60.7 2 <0.0001 2 (R 2GLMM = 0.14) Ovipositing females per census Microhabitat type 6.1 2 0.0474 Species 1.3 1 0.2475 Microhabitat type  species 21.2 2 <0.0001 Note: Results from the type III Wald χ analysis of deviance of the fixed effects. Marginal R for GLMMs (including only fixed effects) were calculated following 2 2 Nakagawa et al. (2017) with the version 0.9.2 of the performance package (Lüdecke et al., 2021). range of 5–10 C in the leaf underside, where most eggs most variable daily thermal profiles (Figure 2g), with were laid, Figure 2b). The thermal regimes characterized more positively skewed distributions (i.e., more extreme throughout the monitoring campaign varied between values in the upper side of the thermal distribution, microhabitats in the same line. Daily maximum tempera- Appendix S1: Figure S8a). Thermal heterogeneity at fine tures recorded with the data loggers were, respectively, a spatial scales was also higher in the open microhabitats, mean of 6 and 2 C lower in the semi-closed and which was calculated as the SD of the temperatures semi-open microhabitats than in the open areas, where recorded on the upper- and underside of the leaves in P. rapae oviposits (Figure 2c). A similar pattern was each microhabitat on a daily basis (Figure 2h). found for temperatures recorded at the upper surface of the leaves (Figure 2d), which were measured during the period of highest insolation (10–16 h) and when most Sources of thermal heterogeneity at fine eggs are usually laid (Appendix S1: Figure S4). scales Canopy cover of the closed and semi-closed micro- habitats from the mid-elevation site diminished the mean All the thermal variables measured at the host-plant level daily temperatures of the understory by an average (i.e., excluding air temperatures at the microhabitat scale of 3 C relative to mean macroclimatic conditions recorded with the data loggers) were strongly correlated (Figure 2e). In contrast, the buffering capacity of lowland in pairwise correlations for each site and microhabitat vegetation was much lower, and the largest offsets from type (R 2 = 0.81 ± 0.16; and p = 0.0042 ± 0.0198, which macroclimatic mean temperatures were observed in the summarize the mean ± sd of the R 2 and p values open and semi-open microhabitats (daily mean tempera- obtained in the correlation tests). Nevertheless, there tures 1–2 C higher than meteorological records). Besides were notable differences between leaves, soils, and the the higher air temperatures at the microhabitat scale, air in the open microhabitat of the lowland site, where leaves in the open microhabitats of the lowland site were thermal heterogeneity was higher (Figure 3 and additionally subject to important thermal amplification Appendix S1: Figures S5 and S9). processes that could elevate foliar temperatures up to The soil in the open microhabitat of the lowland site 10 C higher than the air above the host plant (Figure 2f). was warmer and drier than in the other microhabitats Overall, the structure of the open–closed ecotones (Figures 3 and 4a), especially during summer (Figure 4c). largely determined the microclimatic processes that Soil surface reached temperatures higher than 45 C when operated in the different microhabitats, modifying the exposed to full radiation (i.e., soil thermal amplification, thermal exposure of the host plants. Microhabitat and Figure 3c). Going up from the soil, air temperature rap- foliar temperatures could exceed 40 C in the open idly decreased with height following a hyperbolic sine microhabitats (Figure 2b–d), especially in the lowland function (i2 in Figure 3). In this line, basal leaves of the site during summer (Appendix S1: Figure S5), while tem- host plants from open microhabitats were 2–7 C warmer peratures remained below 40 C in semi-closed and than apical leaves (i1 in Figure 3), while basal leaves in semi-open microhabitats of both sites across the seasons semi-open and semi-closed microhabitats could reach (Appendix S1: Figures S6 and S7). inferior temperatures than apical leaves (i.e., soil cooling Temperatures in closed and semi-closed microhabi- effect, Figure 3b). This thermal difference between leaves tats were more constant throughout the day and were could be recorded in spring plants of Lepidium draba, more homogenous across this type of microsite which were >40 cm in height. In contrast, leaves of sum- (Figure 2g,h). Open microhabitats presented instead the mer resprouts (<5 cm height) presented the highest 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ECOLOGICAL MONOGRAPHS 9 of 19 temperatures and the lowest thermal heterogeneity (daily ranges of microhabitat air temperature (i.e., 20–35 C), SD of foliar thermal records, Appendix S1: Figure S5c,d). foliar underside temperatures were 1–3 C cooler than Fine-scale thermal heterogeneity was also detected upper parts of the leaves exposed to direct radiation in between the sides of the same leaves. In intermediate open microhabitats (Figure 3d). However, at higher air (a) (b) 1.0 40 * of ovipositions temperature Oviposition frequency Relative 35 (ºC) 0.5 30 25 0.0 C SC+SO O P. napi P. rapae Microhabitat Butterfly species (c) (d) a a b c a ab b c 40 Temperature 40 Daily Tmax (ºC) (ºC) Foliar thermal conditions 30 20 20 a b c d a a b b C SC SO O C SC SO O Microhabitat Microhabitat (e) (f) 5 Microhabitat thermal conditions a b c d ab a b c (sensor–meteo,ºC) Thermal offset Thermal offset (leaf–air,ºC) 0 10 –5 0 a b c d a a b b C SC SO O C SC SO O Microhabitat Microhabitat (g) (h) a b c d a a a b Spatial thermal a ab c bc (Daily SD, ºC) (Daily SD, ºC) heterogeneity Daily thermal variability 10 4 2 0 a b c c 0 C SC SO O C SC SO O Microhabitat Microhabitat Microhabitat P. napi P. rapae Site Mid-elevation Lowland selected by FIGURE 2 Legend on next page. 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 10 of 19 VIVES-INGLA ET AL. temperatures (i.e., >35 C), thermal differences between Interspecific differences in thermal the foliar upper and underside vanished (Figure 3e). In strategies and mortality these conditions, leaves could be 10 C warmer than the air temperature above the host plant (Figure 3f). The Time to larval death and experimental temperatures were effects of foliar height and foliar side on thermal hetero- associated in semilogarithmic curves for both species geneity were putatively associated with radiative heating (Figure 5a), presenting a steeper slope (25% more nega- and sensible heat fluxes from the soil and with processes tive) for P. napi (species  temperature p = 0.0085 in the of stomatal closure of the leaves. ANCOVA model, which explained a 76% of the variance). The estimated intersection point of the two TDT curves was located at T* = 41.2 C. Above this threshold, Host-plant variability P. rapae exhibited greater survival than P. napi, while the opposite was true for inferior temperatures (Figure 5a). Open–closed ecotones also induced changes in As predicted, the estimates (±SE) of thermal tolerance host-plant traits, paralleling variation in thermal expo- z and CTmax were higher for P. rapae (z = 5.10 ± 0.24 C; sure (Figure 4 and Appendix S1: Figures S6 and S7). CTmax = 53.48 C) than for P. napi (z = 4.10 ± 0.26 C; Host plants in open microhabitats had smaller leaves CTmax = 51.08 C). In the general mixed linear models with lower ratios of water content (i.e., less water per considering more factors than temperature, strong effects mg of foliar dry weight, Figure 4d). Variations in foliar were also found for the larval weight (Table 2), with chlorophyll content and fruit production depended on smaller larvae presenting lower survival times. the shade tolerance of the host plant. Lepidium draba Concretely, an increase of 0.1 g in larval weight would (lowland site) presented higher chlorophyll contents cancel out the effects of an increase of 1 C in temperature and produced more fruits in open microhabitats, while (i.e., the former would provoke an increase of 30%–50% plants in closed and semi-closed microhabitats did not in survival times, while the latter would decrease them reproduce sexually and had thinner leaves, with low by 30%–50%, Appendix S1: Figure S13). chlorophyll contents. On the contrary, chlorophyll con- Interspecific differences in the TDT curves and the tent and fruit production for A. petiolata in the predicted thermal mortality in the field were in agree- mid-elevation site were lowest in the open microhabitat ment with the pattern of microhabitat selection of the (Appendix S1: Figure S10). two species. Daily mortality was higher for the Foliar water and chlorophyll contents decreased in thermosensitive P. napi (Appendix S1: Figure S14) when both host plants (Figure 4e,f, and Appendix S1: thermal stresses were more acute (more intense heat Figure S11) as they senesced after fructification during late challenges, mainly found in open microhabitats). spring and early summer (ordinal days 140–180). Only In contrast, more constant thermal regimes, with less non-flowering first-year rosettes (A. petiolata) and summer extreme but longer thermal heat challenges (i.e., chronic rhizome resprouts (L. draba) remained in midsummer stress), were deadlier for P. rapae. Daily thermal mortal- after senescence (Appendix S1: Figure S12). First-year ity from March to September was usually low (<0.01), rosettes of A. petiolata notably coexisted in June with the although it could reach values around 0.4 in the open reproductive stage of second-year individuals. In contrast, microhabitats during the warmest days (<15% of the there was no temporal overlap between reproductive days, Appendix S1: Figure S15a). L. draba plants and new summer resprouts, leading to a The accumulation of low daily mortalities for a period of scarcity of fresh host plants lasting 2–3 weeks. period similar to larval development time (i.e., 30 days) F I G U R E 2 Interspecific differences in the thermal regime that eggs and larvae are exposed to, both at microhabitat (c, e and g) and foliar (b, d, f and h) scales. (a) Relative distribution of ovipositions for each species across the open–closed ecotones. (b) Foliar temperature during oviposition at the underside. (c) Daily maximum temperatures recorded with the data loggers. (d) Foliar temperatures at the upper side during host-plant monitoring. (e) Thermal offset calculated as the difference of the mean daily temperatures between the microhabitat and the macroclimate. (f) Foliar thermal offset calculated as the instantaneous difference in foliar and air temperatures above the host plant at 1 m height. (g) Daily temporal variability (standard deviation, SD) of the temperatures recorded at the microhabitat scale. (h) Thermal heterogeneity of foliar temperatures in the same microhabitat and time. Different letters indicate the microhabitats with p < 0.05 in pairwise Tukey HSD tests of the response variable for each site (light gray for the mid-elevation site, and dark gray for the lowland site). The lower and upper hinges of the box represent the 1st and the 3rd quartiles respectively (Q1, Q3); its inner line, the median; and the length of the box, the IQR = Q3  Q1. The lower whisker represents the smallest value ≥Q1  1.5  IQR; and the upper whisker, the biggest value ≤Q3 + 1.5  IQR. Outliers are not shown. Colored areas in panels C–H indicate the microhabitats selected by each species (i.e., blue: P. napi; orange: P. rapae) to facilitate interspecific comparisons. C, closed; SC, semi-closed; SO, semi-open; and O, open microhabitat. 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ECOLOGICAL MONOGRAPHS 11 of 19 Microhabitat of P. napi offspring Microhabitat of P. rapae offspring (a) Microhabitat and host-plant selection More homogeneous More heterogeneous thermal environment, thermal environment, with buffered extremes with amplified extremes AIR 3–10ºC TEMPERATURES FOLIAR 2–10ºC TEMPERATURES 0ºC 1–3ºC 0–2ºC 2–7ºC SOIL (b) (d) (e) (f) (c) TEMPERATURES SOIL HUMIDITY DROUGHT URE EM PERAT SOIL T Temperature Low High (b) Soil cooling effect Soil thermal amplification (c) Height from soil (cm) 100 5 (i1) (i2) R = 0.86 2 Basal–apical ∆T = 0–2ºC ∆T = 2–7ºC leaves (ºC) P < 0.0001 50 0 R2 = 0.12 P = 0.05 0 30 45 Air temperature (ºC) Underside cooling effect (e) 10 R2 = 0.15 (i3) leaf sides (ºC) Upper–under P < 0.0001 ∆T = 1–3ºC Heat dissipation (d) Canopy thermal 0 Net homogeneization Evapotranspirative radiation cooling 25 35 ∆T = 0ºC Air Tmax (ºC) C 0º ∆T = 0ºC R2 = 0.49 Leaf – Air (ºC) (i4) P < 0.0001 10 0ºC Reduced 2–1 evapotranspiration 0 Foliar thermal amplification (f) FIGURE 3 Legend on next page. 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 12 of 19 VIVES-INGLA ET AL. Both sites Mid-elevation Lowland (a) (b) (c) b c bc a 60 60 ab c b a 60 humidity humidity humidity 40 40 40 Soil (%) Soil Soil (%) (%) 20 20 20 0 0 C SC SO O 60 120 180 240 60 120 180 240 300 Microhabitat Ordinal day Ordinal day (d) (e) (f) c b b a 10 b b a a 10 Foliar water Foliar water Foliar water (g g–1 DW) (g g–1 DW) (g g–1 DW) content content content 5 5 5 0 0 C SC SO O 60 120 180 240 60 120 180 240 300 Microhabitat Ordinal day Ordinal day Mid-elevation Microhabitat O SO SC Site Lowland selected by P. rapae P. napi F I G U R E 4 Spatial and seasonal variation in the soil humidity (a–c) and foliar water content (d–f) of the host plants measured during the monitoring campaign. (a and d) Variation across open–closed ecotones. Different letters indicate the microhabitats with p < 0.05 in pairwise Tukey HSD tests of the response variable for each site (light gray for the mid-elevation site, and dark gray for the lowland site). (b, c, e and f) Seasonal cycle of the soil humidity (b and c) and the foliar water content (e and f) measured in different microhabitats types monitored in the mid-elevation (b and e) and lowland (c and f) sites. Open microhabitats (mainly selected by P. rapae) are represented in yellow; and semi-open and semi-closed microhabitats (mainly selected by P. napi), in blue. C, closed; SC, semi-closed; SO, semi-open; and O, open microhabitats. F I G U R E 3 A visual summary of the diverse microclimatic mosaics and processes at fine scales quantified in the lowland site. Blue areas (left) indicate microhabitats preferentially selected by P. napi, and orange areas (right), by P. rapae. Observed thermal differences between the air, host plants, and soil allowed the identification of six processes determining larval thermal exposure: Butterfly oviposition behavior (a); the influence of soil cooling effects (b) and soil thermal amplification on basal leaves (c); canopy thermal homogenization (d); leaf underside cooling by active stomatal conductance (e); and foliar thermal amplification by reduced stomatal conductance (f). Detailed thermal data are reported in the four insets providing evidence for these processes (i1–i4). (a) Air temperatures in the semi-open and semi-closed microhabitats where P. napi oviposits reach about 3–10 C inferior values and are less variable than in open microhabitats (ovipositing microsites of P. rapae). Open microhabitats are more thermally heterogeneous at fine scales, with notable thermal differences between leaves, soil, and the air. (b) Temperatures at the soil surface influence the heat balance of basal leaves, such that basal leaves in the semi-closed microhabitat can be cooler than apical leaves (i1). (c) Soil thermal amplification in the open microhabitat warms the lower air layers and creates steep thermal gradients with height (i2), which can be also detected between basal and apical leaves (i1). (d) Canopy shading in more closed microhabitats cools leaves and air, resulting in similar temperatures in the leaf upper- and undersides (i3 and i4, results in blue). (e) When temperatures in open microhabitats were moderately warm (around 25 C), we observed that stomatal conductance and evapotranspiration cooled foliar undersides 1–3 C relative to the upper sides, which are more directly exposed to radiative heating (i3, orange line). (f) However, in dry and warm conditions (>35 C), high radiative heating from the soil and reduced evapotranspiration linked to leaf stomatal closure can bring leaves 10 C warmer than the air (i4, orange). Black lines in the soil are added for illustrative purposes and represent the patterns of variation of soil temperature and humidity across the ecotone that we observed in Figure 4a and Appendix S1: Figure S11b. More detailed results of thermal heterogeneity at fine scales can be found in Appendix S1: Figures S5 and S9. 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ECOLOGICAL MONOGRAPHS 13 of 19 (a) Species P. napi P. rapae 420 (min, log scale) Survival time 180 60 10 T* 40 42 44 Temperature (°C) (b) Semi-closed and semi-open microhabitats Open microhabitats 0.20 0.7 Mean thermal 0.15 mortality 0.6 0.10 0.5 0.05 0.4 0.3 (c) 2.0 2.0 (P. napi/P. rapae) thermal mortality Ratio of mean 1.5 1.5 1.0 1.0 0.5 0.5 SO C SO O O −O −O O O O − d −S − −S −S e e Ld− Ld− Ld− − Me −L − Me −L d −L d −M −M 8− 9− 0− s1 s2 s3 s4 s5 s6 s7 s s s1 Sensor − site − microhabitat F I G U R E 5 (a) Thermal death time (TDT) curves for P. napi and P. rapae. The dashed lines indicate the intersection point between the two TDT curves, representing the temperature at which both species show equal survival times. (b) Interspecific differences of mean thermal mortality during development for the periods where 30-day cumulative mortality was equal to or higher than 0.01. Colored areas represent the density functions of the mean thermal mortality in the 100 bootstrap replicates for each species and each data logger (each temporal series of microclimatic temperatures). The lower and upper limits of the boxes represent the 1st and the 3rd quartiles of mean thermal mortality and, their inner line, the median. (c) Interspecific ratio (P. napi/P. rapae) of mean thermal mortality during development observed in each microhabitat. For both species, estimated mortality was lower in their preferred microhabitats: Pieris napi showed inferior mortalities in semi-open and semi-closed microhabitats (left panels) but higher in open microhabitats (right panels). Microhabitats with cumulative mortality inferior to 0.01 during the whole period (March–September) are not shown. Me, mid-elevation; Ld, lowland sites. C, closed; SC, semi-closed; SO, semi-open; and O, open microhabitats. 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 14 of 19 VIVES-INGLA ET AL. TABLE 2 General linear mixed model of the time to larval death for each species. Model Fixed effect χ2 df p  Pieris napi (R 2 GLMM = 0.75) Temperature ( C) 194.0 1 <0.0001 Larval weight (g) 5.7 1 0.0166 Site 1.6 1 0.2051 Pieris rapae (R 2GLMM = 0.8) Temperature ( C) 476.1 1 <0.0001 Larval weight (g) 9.4 1 0.0022 Site 1.9 1 0.1691 Note: Results from the type III Wald χ analysis of deviance of the fixed effects. Marginal R for GLMMs (including only fixed effects) were calculated following 2 2 Nakagawa et al. (2017) with the version 0.9.2 of the performance package (Lüdecke et al., 2021). could exert considerable thermal pressure on natural Kingsolver, 1979; Potter et al., 2009; Woods et al., 2022). populations (cumulative mortality >0.8 during the sum- For the first time, here we predicted the mortality derived mer in the open microhabitats). However, the estimates from microclimatic temperatures considering both the indicated that thermal mortality would be importantly intensities and the duration of the thermal exposures and reduced if larvae conducted thermal avoidance behaviors their cumulative effects (Rezende et al., 2020b). (Appendix S1: Figure S15b). Mean thermal mortality dur- The two butterflies selected the same host-plant species ing the development was around 0–0.2 in semi-open and from different microhabitats to oviposit: P. napi laid most semi-closed microhabitats (Figure 5b), and an average of of their eggs in semi-open and semi-closed microhabitats, 15% lower for P. napi than for P. rapae (Figure 5c). The and P. rapae, in open ones (Table 1 and Figure 2). These opposite pattern was found in the open microhabitats, results indicate that habitat choice preceded host-plant where mean thermal mortality during development was selection, as has long been proposed (Courtney, 1986; 15% higher for P. napi than for P. rapae and ranged Dennis, 2010; Porter, 1992) and as has been found for this between 0.4 and 0.8. Thus, by laying its eggs on pair of species (Friberg & Wiklund, 2019; Ohsaki, 1982; semi-open and semi-closed microhabitats, P. napi eludes Ohsaki & Sato, 1999). Ovipositing females generally use a high thermal mortality. In contrast, P. rapae selects the multiple cues following a spatially-structured and hierar- microhabitats with the most intense heat challenges, chical process, from coarse to finer scales. Determining exposing its offspring to deadly thermal stresses that it which specific cues influenced oviposition decisions is can withstand better than P. napi. beyond the scope of our study, but both microclimatic and host-plant factors could have had a role. Ohsaki (1982) associated P. napi and P. rapae oviposition decisions with DISC USS I ON the different light conditions of the microhabitat, and Forsberg (1987) suggested that P. napi in Sweden actively Our work quantified the microclimatic and host-plant oviposited in small plants to favor higher microclimatic conditions of the microhabitats where two butterfly spe- temperatures. Other studies related oviposition decisions cies lay their eggs and assessed their potential impacts on of P. rapae with host-plant qualities. Visual stimuli the offspring. We conducted the study with two model (mainly the color or the greenness of the plant) have a key species from an extensively-studied family of butterflies role when females are searching for a host plant in the ecotones between forested and open habitats of (Myers, 1985; Tsuji & Coe, 2014), although olfactory cues two protected areas. We combined field censuses of but- likely influence pre-alightment decisions too (Renwick & terfly behavior, detailed measurements of the microcli- Radke, 1988). After landing on host plants, the chemical mate in the selected microhabitats, and an experimental and nutrient status of the plants were found to be decisive characterization of the thermal tolerance of the larvae, to (Hern et al., 1996). For example, leaves with higher water computationally predict the thermal mortality in the and nitrogen content, and with high transpiration rates, field. Many studies stress the need to consider tempera- are more frequently accepted to oviposit (Myers, 1985; tures at the microclimatic scale to better understand the Wolfson, 1980). In our study, all of these factors varied thermal ecology of the species and their responses to between the selected microhabitats (Figures 2 and 4). global warming (Bramer et al., 2018; Pincebourde & Microhabitat preferences observed in the field Woods, 2020; Woods et al., 2015). As far as we know, matched interspecific differences in the TDT curves mortality associated with microclimatic variation has determined in the laboratory and the predicted thermal been very occasionally assessed (Kaiser et al., 2016; mortality in the field (Figure 5). P. napi presented a 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ECOLOGICAL MONOGRAPHS 15 of 19 more thermosensitive strategy, with lower survivals where soil cooling and canopy shading homogenize tem- under acute stresses, but higher under longer peratures (Figure 3b,d). Our predictions of thermal mor- subextreme challenges. Accordingly, its mean thermal tality were based on the thermal series extracted from the mortality during development was approximately 15% data loggers, which represent air temperatures at the higher than P. rapae in open microhabitats, where heat microhabitat level. But larvae might be more dependent challenges are more extreme, but 15% lower than on temperatures at the foliar level, especially at their ini- P. rapae in their selected microhabitats. Petersen (1954) tial stages (Pincebourde et al., 2021; Woods, 2013). The already proposed that, as a result of interspecific compe- high fine-scale thermal heterogeneity in open microhabi- tition, P. napi would have specialized in developing in tats could potentially expose larvae of P. rapae to more shaded environments, where thermal regimes are cooler acute thermal stresses than those we predicted, but could and host plants have thinner leaves than in the dry, also offer more opportunities for behavioral thermoregu- open habitats selected by P. rapae. But later research lation (by moving to the foliar underside [Figure 3e], or a suggested that competition between Pierids is unlikely cooler leaf or plant [Figure 3c]). Although thermal het- to be driving their habitat segregation (Courtney, 1986), erogeneity was lower in summer in the lowland, some and that escape from parasitism is a more likely driver leaves presented temperatures that could reduce thermal in the case of P. napi and P. rapae (Ohsaki & mortality. For example, if larvae avoided thermal expo- Sato, 1999). Regardless of the ultimate cause driving this sures >40 C, cumulative mortality during development habitat differentiation, specialization in fresh and would be importantly reduced (Appendix S1: Figures S5 humid habitats or dry and hot habitats usually comes and S15). Thermal avoidance behaviors have been with different costs and benefits. For example, studies reported for the larvae of P. napi (Carnicer et al., 2019) with other butterflies have shown that specializing in and P. rapae (Kingsolver & Gomulkiewicz, 2003). hot environments can select higher fecundities and Microclimatic exposure during the development of adult survival at the expense of larval survival insects will simultaneously affect other physiological and (Karlsson & Wiklund, 2005). In this line, our results also demographical rates besides larval survival (Braschler predicted that larval survival in open microhabitats is et al., 2021; Diamond et al., 2013; Kingsolver et al., 2011; lower (Figure 5), and previous studies found that von Schmalensee et al., 2021), which can influence the P. rapae laid more eggs but smaller ones than P. napi patterns of larval mortality described here. For example, (Ohsaki, 1982). different microclimatic regimes can lead to shorter or In this study, we also identified six processes longer larval development times due to nonlinear (Figure 3a–f) modulating the thermal exposure of the two changes in development rates with temperature (thermal species at very fine scales, and hence their thermal mor- performance curves; Greiser et al., 2022; von Schmalensee tality. The first one is microhabitat selection by adult but- et al., 2021; Appendix S1: Figure S16). These variations in terflies, which strongly determines the thermal exposure the development time would in turn modulate the time of the larvae. Thermal differences between open and of exposure to microclimatic lethal stresses. Moreover, intermediately-covered microhabitats were in the range our results indicated that smaller larvae were more sensi- of 3–10 C for diurnal temperatures at the foliar and air tive to heat challenges (Table 2 and Appendix S1: levels (Figures 2 and 3a). These results are in agreement Figure S13). TDT curves in this study were built with all with the thermal differences between more open and the larvae, so predictions of thermal mortality assumed a closed habitats of herbivorous insects reported in other constant weight during the whole development. studies (Ashton et al., 2009; Friberg et al., 2008; Merckx However, as larvae grow, they would become more toler- et al., 2015; Ohsaki, 1982; Suggitt et al., 2012). Daily tem- ant, suggesting that plastic and evolutionary changes in peratures in open microhabitats also presented higher growth rates and thermal performance curves could play variation, with more extreme temperatures. Then, heat a key role in larval thermal mortality. Further studies and water fluxes at fine scales create very different micro- should consider how growth dynamically modifies ther- climatic mosaics between the selected microhabitats mal tolerance landscapes during development, and the (Figure 3b–f). In open microhabitats, processes like soil parallel effects of microclimatic regimes on growth, and foliar thermal amplification can raise foliar tempera- development, and mortality. tures 10 C relative to the air (Figure 3c,f; Carnicer The differences between the microclimatic regimes of et al., 2021; Pincebourde et al., 2021; Woods et al., 2022), the microhabitats where the two species lay their eggs led or evapotranspiration can cool the foliar underside 3 C to contrasting patterns of thermal mortality of the off- relative to the upper side (Figure 3e). These processes cre- spring. Despite being more thermosensitive, the offspring ate a more heterogeneous thermal mosaic than that of P. napi were predicted to present higher survivals than found in semi-open and semi-closed microhabitats, P. rapae (Figure 5) due to the thermal buffering provided 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 16 of 19 VIVES-INGLA ET AL. by the microhabitats with intermediate cover (Figure 2). Melero et al., 2016; Stefanescu et al., 2011; Ubach Our results feed the growing literature on the key role et al., 2021). Predictive models of ecological responses to that vegetation cover may play in buffering thermal climate change should incorporate information on cli- macroclimatic stresses (Carnicer et al., 2019; De Frenne matic exposure at relevant scales and capture the key et al., 2021; Zellweger et al., 2020). However, many stud- processes shaping organisms’ sensitivity and perfor- ies on butterflies indicate that other species-specific mance in the dynamic thermal conditions they experi- requirements, such as host plant condition and availabil- ence in nature. ity, should be met in these microhabitats to successfully buffer macroclimatic impacts (Ashton et al., 2009; ACKNOWLEDGMENTS Bennett et al., 2015; Carnicer et al., 2019; Kaiser Carlos L opez, Sofía Cortizas, Joaquim de Gispert, et al., 2016; Nieto-Sanchez et al., 2015; Stefanescu Katarzyna Bartnik, and Agnieszka Juszczak contributed et al., 2011; Suggitt et al., 2011, 2012). For example, in to the field and experimental work. Roger Vila summer in the lowland site, L. draba plants reduced provided insightful comments that improved the manu- their foliar water and chlorophyll contents until their script. Emili Bassols provided support with permission complete senescence (Figure 4 and Appendix S1: management, scientific advice and key assistance Figure S11), which would likely interact with the micro- during fieldwork. Parc Natural de la Zona Volcànica climatic impacts on larval development (Clissold & de la Garrotxa and Parc Natural Aiguamolls de Simpson, 2015). We suggest that the different adaptive l’Empordà provided logistic support. This research was strategies of the species and their different oviposition supported by the Spanish Ministry of Science and behaviors pose different challenges to the populations. Innovation through a doctoral grant (FPU17/05869); the With the current microclimatic conditions, P. napi may Spanish Ministry of Economic Affairs and Digital not be as threatened by thermal exposure as P. rapae. Transformation (PID2020-117636GB-C21; CGL2016-78093-R However, the decay and disappearance of L. draba during and CGL2013-48074-P); the Catalan Government the summer in the lowlands might have more negative (SGR-2017-1005); and by Nederlandse Organisatie voor impacts on P. napi because they rely on fewer host-plant Wetenschappelijk Onderzoek (863.11.021). species than P. rapae, which are more likely to find other host plants in a better condition during this period CONFLICT OF INTEREST (Carnicer et al., 2019; Vives-Ingla et al., 2020). The authors declare no conflicts of interest. Our study assessed both the thermal exposure and the thermal sensitivity of two species to predict the associated DA TA AVAI LA BI LI TY S T ATE ME NT mortality using a dynamic model calibrated with physio- Data and code (Vives-Ingla et al., 2022) is available on logical information from the experiments and simulated Zenodo at https://doi.org/10.5281/zenodo.7358091. with the microclimatic regimes recorded in the field. Previous assessments of the vulnerability of organisms to ORCID global warming usually compared their experimental Maria Vives-Ingla https://orcid.org/0000-0003-4887- upper thermal limits with maximum temperatures 8392 recorded in the field (Duffy et al., 2015; Pincebourde & Javier Sala-Garcia https://orcid.org/0000-0003-4545- Casas, 2019; Sunday et al., 2014; Woods et al., 2022). 2580 However, this approach overlooks the time-dependent Constantí Stefanescu https://orcid.org/0000-0001-8952- effects of thermal stress and the cumulative nature of heat 7869 injury (Jørgensen et al., 2021; Rezende et al., 2020b), Josep Peñuelas https://orcid.org/0000-0002-7215-0150 which may underestimate the vulnerability of organisms Jofre Carnicer https://orcid.org/0000-0001-7454-8296 to global warming (Huey & Kearney, 2020). We predicted larval thermal mortalities applying the tolerance landscape RE FER EN CES framework, which allowed us to adopt a more realistic Ashton, S., D. Gutiérrez, and R. J. Wilson. 2009. “Effects of approach by considering their thermal sensitivity to both Temperature and Elevation on Habitat Use by a Rare the duration and the intensity of microclimatic thermal Mountain Butterfly: Implications for Species Responses to exposures (Rezende et al., 2014, 2020b). This new frame- Climate Change.” Ecological Entomology 34: 437–46. Awmack, C. S., and S. R. Leather. 2002. “Host Plant Quality and work can offer new insights into the role of climate Fecundity in Herbivorous Insects.” Annual Review of change in the declines reported for many insect species Entomology 47: 817–44. around the globe (Didham et al., 2020; Wagner Bates, D., M. Mächler, B. Bolker, and S. Walker. 2015. “Fitting et al., 2021, and all the references therein), also in our Linear Mixed-Effects Models Using lme4.” Journal of study area (Colom et al., 2022; Herrando et al., 2019; Statistical Software 67: 1–48. 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ECOLOGICAL MONOGRAPHS 17 of 19 Bennett, N. L., P. M. Severns, C. Parmesan, and M. C. Singer. 2015. Physiology to Predict the Responses of Ants to Climatic “Geographic Mosaics of Phenology, Host Preference, Adult Warming.” Integrative and Comparative Biology 53: 965–74. Size and Microhabitat Choice Predict Butterfly Resilience to Didham, R. K., Y. Basset, C. M. Collins, S. R. Leather, N. A. Climate Warming.” Oikos 124: 41–53. Littlewood, M. H. M. Menz, J. Müller, et al. 2020. Bramer, I., B. J. Anderson, J. J. Bennie, A. J. Bladon, P. De Frenne, “Interpreting Insect Declines: Seven Challenges and a Way D. L. Hemming, R. A. Hill, et al. 2018. “Advances in Forward.” Insect Conservation and Diversity 13: 103–14. Monitoring and Modelling Climate at Ecologically Relevant Doak, P., P. Kareiva, and J. Kingsolver. 2006. “Fitness Scales.” Advances in Ecological Research 58: 101–61. Consequences of Choosy Oviposition for a Time-Limited Braschler, B., S. L. Chown, and G. A. Duffy. 2021. “Sub-Critical Butterfly.” Ecology 87: 395–408. Limits Are Viable Alternatives to Critical Thermal Limits.” Duffy, G. A., B. W. Coetzee, C. Janion-Scheepers, and S. L. Chown. Journal of Thermal Biology 101: 103106. 2015. “Microclimate-Based Macrophysiology: Implications for Carnicer, J., C. Stefanescu, M. Vives-Ingla, C. Lopez, S. Cortizas, Insects in a Warming World.” Current Opinion in Insect C. W. Wheat, R. Vila, J. Llusià, and J. Peñuelas. 2019. Science 11: 84–9. “Phenotypic Biomarkers of Climatic Impacts on Declining Everiit, J. H., A. J. Richardson, and H. W. Guasman. 1985. “Leaf Insect Populations: A Key Role for Decadal Drought, Thermal Reflectance-Nitrogen-Chlorophyll Relations in Buffelgrass.” Buffering and Amplification Effects and Host Plant Photogrametria Engineering and Remote Sensing 51: 463–6. Dynamics.” Journal of Animal Ecology 88: 376–91. Forsberg, J. 1987. “Size Discrimination among Conspecific Carnicer, J., M. Vives-Ingla, L. Blanquer, X. Méndez-Camps, Hostplants in Two Pierid Butterflies; Pieris napi L. and Pontia C. Rosell, S. Sabaté, E. Gutiérrez, T. Sauras, J. Peñuelas, and daplidice L.” Oecologia 72: 52–7. A. Barbeta. 2021. “Forest Resilience to Global Warming Is Francis, A., and S. I. Warwick. 2008. “The Biology of Canadian Strongly Modulated by Local-Scale Topographic, Weeds. 3. Lepidium draba L., L. chalepense L., L. appelianum Microclimatic and Biotic Conditions.” Journal of Ecology 109: Al-Shehbaz.” Canadian Journal of Plant Science 88: 379–401. 3322–39. Franken, O., M. Huizinga, J. Ellers, and M. P. Berg. 2018. “Heated Carnicer, J., C. W. Wheat, M. Vives-Ingla, A. Ubach, Communities: Large Inter- and Intraspecific Variation in Heat C. Domingo-Marim on, S. Nylin, C. Stefanescu, R. Vila, Tolerance across Trophic Levels of a Soil Arthropod C. Wiklund, and J. Peñuelas. 2017. “Evolutionary Responses of Community.” Oecologia 186: 311–22. Invertebrates to Global Climate Change: The Role of Friberg, M., M. Olofsson, D. Berger, B. Karlsson, and C. Wiklund. Life-History Trade-Offs and Multidecadal Climate Shifts.” In 2008. “Habitat Choice Precedes Host Plant Choice - Niche Global Climate Change and Terrestrial Invertebrates, edited by Separation in a Species Pair of a Generalist and a Specialist S. N. Johnson and H. Jones, 319–48. Chichester: John Wiley & Butterfly.” Oikos 117: 1337–44. Sons, Ltd. Friberg, M., and C. Wiklund. 2019. “Host Preference Variation can- Cavers, P. B., M. I. Heagy, and R. F. Kokron. 1979. “The Biology of not Explain Microhabitat Differentiation among Sympatric Canadian Weeds. 35. Alliaria petiolata (M. Bieb.) Cavara and Pieris napi and Pieris rapae Butterflies.” Ecological Entomology Grande.” Canadian Journal of Plant Science 59: 217–29. 44: 571–6. Clissold, F. J., and S. J. Simpson. 2015. “Temperature, Food Quality Geiger, R. 1950. The Climate near the Ground, 2nd ed. Cambrigde: and Life History Traits of Herbivorous Insects.” Current Harvard University Press. Opinion in Insect Science 11: 63–70. Gibbs, M., and H. Van Dyck. 2009. “Reproductive Plasticity, Colom, P., M. Ninyerola, X. Pons, A. Traveset, and C. Stefanescu. Oviposition Site Selection, and Maternal Effects in 2022. “Phenological Sensitivity and Seasonal Variability Fragmented Landscapes.” Behavioral Ecology and Sociobiology Explain Climate-Driven Trends in Mediterranean Butterflies.” 64: 1–11. Proceedings of the Royal Society B: Biological Sciences 289: Greiser, C., L. von Schmalensee, O. Lindestad, K. Gotthard, 20220251. and P. Lehmann. 2022. “Microclimatic Variation Affects Courtney, S. P. 1986. “The Ecology of Pierid Butterflies: Developmental Phenology, Synchrony and Voltinism in an Dynamics and Interactions.” Advances in Ecological Insect Population.” Functional Ecology 36: 3036–48. Research 15: 51–131. Hern, A., G. Edwards-Jones, and R. G. McKinlay. 1996. “A Review Curran, P. J., J. L. Dungan, and H. L. Gholz. 1990. “Exploring the of the Pre-Oviposition Behaviour of Small Cabbage White Relationship between Reflectance Red Edge and Chlorophyll Butterfly, Pieris rapae (Lepidoptera: Pieridae).” Annals of Content in Slash Pine.” Tree Physiology 7: 33–48. Applied Biology 128: 349–71. de Bolos, O., and J. Vigo. 1990. Flora dels Països Catalans, II. First. Herrando, S., N. Titeux, L. Brotons, M. Anton, A. Ubach, D. Villero, Barcino, Barcelona: Pòrtic. E. García-Barros, M. L. Munguira, C. Godinho, and C. De Frenne, P., J. Lenoir, M. Luoto, B. R. Scheffers, F. Zellweger, Stefanescu. 2019. “Contrasting Impacts of Precipitation on J. Aalto, M. B. Ashcroft, et al. 2021. “Forest Microclimates and Mediterranean Birds and Butterflies.” Scientific Reports 9: 5680. Climate Change: Importance, Drivers and Future Research Huey, R. B., and M. R. Kearney. 2020. “Dynamics of Death by Agenda.” Global Change Biology 27: 2279–97. Heat.” Science 369: 1163. Dennis, R. L. H. 2010. A Resource-Based Habitat View for Jørgensen, L. B., H. Malte, M. Ørsted, N. A. Klahn, and Conservation: Butterflies in the British Landscape. West Sussex: J. Overgaard. 2021. “A Unifying Model to Estimate Thermal John Wiley & Sons, Ltd. Tolerance Limits in Ectotherms across Static, Dynamic and Diamond, S. E., C. A. Penick, S. L. Pelini, A. M. Ellison, N. J. Fluctuating Exposures to Thermal Stress.” Scientific Reports Gotelli, N. J. Sanders, and R. R. Dunn. 2013. “Using 11: 12840. 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 18 of 19 VIVES-INGLA ET AL. Kaiser, A., T. Merckx, and H. Van Dyck. 2016. “The Urban Heat Area III. Difference in the Annual Generation Numbers in Island and its Spatial Scale Dependent Impact on Survival and Relation to Habitat Selection by Adults.” Researches on Development in Butterflies of Different Thermal Sensitivity.” Population Ecology 24: 193–210. Ecology and Evolution 6: 4129–40. Ohsaki, N., and Y. Sato. 1999. “The Role of Parasitoids in Karlsson, B., and C. Wiklund. 2005. “Butterfly Life History and Evolution of Habitat and Larval Food Plant Preference by Temperature Adaptations; Dry Open Habitats Select for Three Pieris Butterflies.” Researches on Population Ecology Increased Fecundity and Longevity.” Journal of Animal 41: 107–19. Ecology 74: 99–104. Petersen, B. 1954. “Egg-Laying and Habitat Selection in some Pieris Kaspari, M., N. A. Clay, J. Lucas, S. P. Yanoviak, and A. Kay. 2015. Species (Lep.).” Entomologisk tidskrift 75: 194–203. “Thermal Adaptation Generates a Diversity of Thermal Limits Pincebourde, S., and J. Casas. 2019. “Narrow Safety Margin in the in a Rainforest Ant Community.” Global Change Biology 21: Phyllosphere during Thermal Extremes.” Proceedings of the 1092–102. National Academy of Sciences 116: 5588–96. Kingsolver, J. G. 1979. “Thermal and Hydric Aspects of Pincebourde, S., M. E. Dillon, and H. A. Woods. 2021. “Body Size Environmental Heterogeneity in the Pitcher Plant Mosquito.” Determines the Thermal Coupling between Insects and Plant Ecological Monographs 49: 357–76. Surfaces.” Functional Ecology 35: 1424–36. Kingsolver, J. G. 2000. “Feeding, Growth, and the Thermal Pincebourde, S., C. C. Murdock, M. Vickers, and M. W. Sears. Environment of Cabbage White Caterpillars, Pieris rapae L.” 2016. “Fine-Scale Microclimatic Variation Can Shape the Physiological and Biochemical Zoology 73: 621–8. Responses of Organisms to Global Change in both Natural Kingsolver, J. G., and R. Gomulkiewicz. 2003. “Environmental and Urban Environments.” Integrative and Comparative Variation and Selection on Performance Curves.” Integrative Biology 56: 45–61. and Comparative Biology 43: 470–7. Pincebourde, S., and H. A. Woods. 2012. “Climate Uncertainty on Kingsolver, J. G., H. A. Woods, L. B. Buckley, K. A. Potter, H. J. Leaf Surfaces: The Biophysics of Leaf Microclimates and their MacLean, and J. K. Higgins. 2011. “Complex Life Cycles and Consequences for Leaf-Dwelling Organisms.” Functional the Responses of Insects to Climate Change.” Integrative and Ecology 26: 844–53. Comparative Biology 51: 719–32. Pincebourde, S., and H. A. Woods. 2020. “There Is Plenty of Room Lenth, R. 2020. “emmeans: Estimated Marginal Means, aka at the Bottom: Microclimates Drive Insect Vulnerability to Least-Squares Means.” https://CRAN.R-project.org/package= Climate Change.” Current Opinion in Insect Science 41: emmeans. 63–70. Lüdecke, D., M. S. Ben-Shachar, I. Patil, P. Waggoner, and Pollard, E., and T. J. Yates. 1993. Monitoring Butterflies for Ecology D. Makowski. 2021. “performance: An R Package for and Conservation: The British Butterfly Monitoring Scheme, 1st Assessment, Comparison and Testing of Statistical Models.” ed. London: Chapman & Hall. Journal of Open Source Software 6: 3139. Poorter, H., Ü. Niinemets, N. Ntagkas, A. Siebenkäs, M. Mäenpää, Melero, Y., C. Stefanescu, and J. Pino. 2016. “General Declines in S. Matsubara, and T. L. Pons. 2019. “A Meta-Analysis of Plant Mediterranean Butterflies over the Last Two Decades Are Responses to Light Intensity for 70 Traits Ranging from Modulated by Species Traits.” Biological Conservation 201: Molecules to Whole Plant Performance.” New Phytologist 223: 336–42. 1073–105. Merckx, T., M. Serruys, and H. Van Dyck. 2015. “Anthropogenic Porter, K. 1992. “Eggs and Egg-laying.” In Ecology of Butterflies in Host Plant Expansion Leads a Nettle-Feeding Butterfly out Britain, edited by R. L. H. Dennis, 46–72. Oxford: Oxford of the Forest: Consequences for Larval Survival and University Press, Inc. Developmental Plasticity in Adult Morphology.” Evolutionary Potter, K. A., G. Davidowitz, and H. A. Woods. 2009. “Insect Eggs Applications 8: 363–72. Protected from High Temperatures by Limited Myers, J. H. 1985. “Effect of Physiological Condition of the Host Homeothermy of Plant Leaves.” Journal of Experimental Plant on the Ovipositional Choice of the Cabbage White Biology 212: 3448–54. Butterfly, Pieris rapae.” Journal of Animal Ecology 54: 193–204. R Core Team. 2019. R: A Language and Environment for Statistical Nakagawa, S., P. C. D. Johnson, and H. Schielzeth. 2017. “The Computing. Vienna, Austria: R Foundation for Statistical Coefficient of Determination R 2 and Intra-Class Correlation Computing. Coefficient from Generalized Linear Mixed-Effects Models Renwick, J. A. A., and C. D. Radke. 1988. “Sensory Cues in Host Revisited and Expanded.” Journal of the Royal Society Interface Selection for Oviposition by the Cabbage Butterfly, Pieris 14: 20170213. rapae.” Journal of Insect Physiology 34: 251–7. Nieto-Sanchez, S., D. Gutiérrez, and R. J. Wilson. 2015. “Long-Term Rezende, E., F. Bozinovic, A. Szilagyi, and M. Santos. 2020a. Change and Spatial Variation in Butterfly Communities over “Dataset and Scripts from: Predicting Temperature Mortality an Elevational Gradient: Driven by Climate, Buffered by and Selection in Natural Drosophila Populations.” https://doi. Habitat.” Diversity and Distributions 21: 950–61. org/10.5061/dryad.stqjq2c1r. Ohsaki, N. 1979. “Comparative Population Studies of Three Pieris Rezende, E. L., F. Bozinovic, A. Szilagyi, and M. Santos. 2020b. Butterflies, P. rapae, P. melete and P. napi, Living in the Same “Predicting Temperature Mortality and Selection in Natural Area: I. Ecological Requirements for Habitat Resources in the Drosophila Populations.” Science 369: 1242–5. Adults.” Researches on Population Ecology 20: 278–96. Rezende, E. L., L. E. Castañeda, and M. Santos. 2014. “Tolerance Ohsaki, N. 1982. “Comparative Population Studies of Three Pieris Landscapes in Thermal Ecology.” Functional Ecology 28: Butterflies, P. rapae, P. melete and P. napi, Living in the Same 799–809. 15577015, 2023, 2, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecm.1561 by JOSEP Penuelas - Readcube (Labtiva Inc.) , Wiley Online Library on [03/05/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ECOLOGICAL MONOGRAPHS 19 of 19 Scriber, J. M., and F. Slansky. 1981. “The Nutritional Ecology of von Schmalensee, L., K. H. Gunnarsd ottir, J. Näslund, K. Gotthard, Immature Insects.” Annual Review of Entomology 26: 183–211. and P. Lehmann. 2021. “Thermal Performance under Stefanescu, C. 2000. “El Butterfly Monitoring Scheme en Constant Temperatures Can Accurately Predict Insect Catalunya: Los primeros cinco años.” Treballs de la Societat Development Times across Naturally Variable Microclimates.” Catalana de Lepidopterologia 15: 5–48. Ecology Letters 24: 1633–45. Stefanescu, C., J. Peñuelas, J. Sardans, and I. Filella. 2006. “Females Wagner, D. L., E. M. Grames, M. L. Forister, M. R. Berenbaum, and of the Specialist Butterfly Euphydryas aurinia (Lepidoptera: D. Stopak. 2021. “Insect Decline in the Anthropocene: Death Nymphalinae: Melitaeini) Select the Greenest Leaves of by a Thousand Cuts.” Proceedings of the National Academy of Lonicera implexa (Caprifoliaceae) for Oviposition.” European Sciences 118: e2023989118. Journal of Entomology 103: 569–74. Wolfson, J. L. 1980. “Oviposition Response of Pieris rapae to Stefanescu, C., I. Torre, J. Jubany, and F. Paramo. 2011. “Recent Environmentally Induced Variation in Brassica nigra.” Trends in Butterfly Populations from North-East Spain and Entomologia Experimentalis et Applicata 27: 223–32. Andorra in the Light of Habitat and Climate Change.” Journal Woods, H. A. 2010. “Water Loss and Gas Exchange by Eggs of of Insect Conservation 15: 83–93. Manduca sexta: Trading off Costs and Benefits.” Journal of Stoutjesdijk, P., and J. J. Barkman. 2014. Microclimate, Vegetation & Insect Physiology 56: 480–7. Fauna. Zeist: KNNV Publishing. Woods, H. A. 2013. “Ontogenetic Changes in the Body Temperature Suggitt, A. J., P. K. Gillingham, J. K. Hill, B. Huntley, W. E. Kunin, of an Insect Herbivore.” Functional Ecology 27: 1322–31. D. B. Roy, and C. D. Thomas. 2011. “Habitat Microclimates Woods, H. A., M. E. Dillon, and S. Pincebourde. 2015. “The Roles of Drive Fine-Scale Variation in Extreme Temperatures.” Oikos Microclimatic Diversity and of Behavior in Mediating the 120: 1–8. Responses of Ectotherms to Climate Change.” Journal of Suggitt, A. J., C. Stefanescu, F. Paramo, T. H. Oliver, B. J. Thermal Biology 54: 86–97. Anderson, J. K. Hill, D. B. Roy, T. M. Brereton, and C. D. Woods, H. A., G. Legault, J. G. Kingsolver, S. Pincebourde, A. A. Thomas. 2012. “Habitat Associations of Species Show Shah, and B. G. Larkin. 2022. “Climate-Driven Thermal Consistent but Weak Responses to Climate.” Biology Letters 8: Opportunities and Risks for Leaf Miners in Aspen Canopies.” 590–3. Ecological Monographs 92: e1544. Sunday, J. M., A. E. Bates, M. R. Kearney, R. K. Colwell, N. K. Zellweger, F., P. De Frenne, J. Lenoir, P. Vangansbeke, Dulvy, J. T. Longino, and R. B. Huey. 2014. “Thermal-Safety K. Verheyen, M. Bernhardt-Römermann, L. Baeten, et al. Margins and the Necessity of Thermoregulatory Behavior 2020. “Forest Microclimate Dynamics Drive Plant Responses across Latitude and Elevation.” Proceedings of the National to Warming.” Science 368: 772–5. Academy of Sciences 111: 5610–5. Zuur, A. F., E. N. Ieno, N. J. Walker, A. A. Saveliev, and G. M. Tsuji, J., and L. Coe. 2014. “Effects of Foliage Color on the Smith. 2009. Mixed Effects Models and Extensions in Ecology Landing Response of Pieris rapae (Lepidoptera: Pieridae).” with R. New York: Springer Science+Business Media. Environmental Entomology 43: 989–94. Ubach, A., F. Paramo, and C. Stefanescu. 2021. “Heterogeneidad En Las Respuestas Demograficas Asociadas Al Gradiente SU PP O R TI N G I N F O RMA TI O N Altitudinal: El Caso De Las Mariposas En El Noreste Ibérico.” Additional supporting information can be found online in Ecosistemas 30: 2148. the Supporting Information section at the end of this article. Vives-Ingla, M., J. Sala-Garcia, C. Stefanescu, A. Casad o-Tortosa, M. Garcia, J. Peñuelas, and J. Carnicer. 2022. “mvives-in- gla/ecotones: Interspecific Differences in Microhabitat Use How to cite this article: Vives-Ingla, Maria, Expose Insects to Contrasting Thermal Mortality (v1.0).” Javier Sala-Garcia, Constantí Stefanescu, Zenodo. https://doi.org/10.5281/zenodo.7358091. Armand Casado -Tortosa, Meritxell Garcia, Vives-Ingla, M., C. Stefanescu, J. Sala-Garcia, and J. Carnicer. Josep Peñuelas, and Jofre Carnicer. 2023. 2020. “Plastic and Phenological Variation of Host Plants “Interspecific Differences in Microhabitat Use Mediates Local Responses of the Butterfly Pieris napi to Expose Insects to Contrasting Thermal Mortality.” Drought in the Mediterranean Basin.” In Butterfly and Moths in l’Empordà and their Response to Global Change, Ecological Monographs 93(2): e1561. https://doi. edited by C. Stefanescu and T. Lafranchis, 113–29. Torroella org/10.1002/ecm.1561 de Montgrí: Recerca i Territori.