Microclimate buffering and thermal tolerance across elevations in a tropical butterfly

Land use and climate change are forcing organisms in the Anthropocene to move, adapt or die (Dirzo et al., 2014). But moving in an increasingly fragmented landscape or adapting to an ever-changing climate might be difficult for organisms usually exposed to a narrow range of environmental conditions. Organisms restricted to stable climates or with limited dispersal abilities have been predicted to be at particular risk of extinction (Bestion et al., 2015; Kingsolver et al., 2013). Despite tropical ectotherms making up half of the animal species described, our knowledge of their potential to cope with high temperatures in natural settings is limited, especially along elevational clines (García-Robledo et al., 2016; Sheldon, 2019). We therefore need a better understanding of the ability of ectotherms to cope with temperatures across elevations and of the climate buffering potential of tropical forests (García-Robledo et al., 2016; Sheldon, 2019).

Since the 1960s, the notion that ‘mountain passes are higher in the tropics’ (Janzen, 1967) has inspired generations of ecological and evolutionary research. Janzen's ‘seasonality hypothesis’ predicts that the reduced seasonality in the tropics selects for narrower thermal tolerances than in temperate zones, which would in turn limit their dispersal across elevations (Angilletta, 2009; Nadeau et al., 2017; Sheldon et al., 2018). Subsequent empirical studies have shown that thermal breadth of insects is indeed higher in temperate zones, where seasonality is stronger than in the tropics (Deutsch et al., 2008; Shah et al., 2018). Furthermore, the great level of specialisation of tropical montane species, reflected by high levels of endemism and beta diversity at high altitudes, may highlight further temperature specificity, and therefore susceptibility to the effects of global warming (Polato et al., 2018; Sheldon, 2019). However, in the face of climate and land-use change in lowland habitats, mountains can act as refugia. Some vulnerable lowland organisms are already shifting their ranges upward (Lawler et al., 2013; Morueta-Holme et al., 2015; Scriven et al., 2015), thus it is critical to ascertain the potential of ectotherms to overcome the physiological barriers that mountains pose.

Janzen's hypothesis has often been tested at the macroecological scale and with interpolated data from weather stations, assuming that tropical ectotherms live at ambient air temperature (Pincebourde and Suppo, 2016). However, this ignores the microclimate differences most relevant to organisms inhabiting tropical forests (Potter et al., 2013). Tropical forests are very heterogeneous habitats with a particularly steep vertical climatic gradient, such that the understorey is often more than 2°C cooler than the canopy and spans an 11% difference in relative humidity (Scheffers et al., 2013). This complexity is not fully captured by interpolated datasets often used in ecological modelling (e.g. WorldClim2; Fick and Hijmans, 2017), with mean temperatures in some cases overestimated by 1.5°C (Blonder et al., 2018; Kearney and Porter, 2009; Storlie et al., 2014). Thus, the biological relevance of studies in the tropics using weather station data is limited, as they are positioned specifically to minimise habitat characteristics that can be crucial in determining the thermal tolerance of local organisms (De Frenne and Verheyen, 2016; Jucker et al., 2018; Senior et al., 2017). These biases could become even more pronounced at higher elevations, where weather stations are very sparse in the tropics (Fick and Hijmans, 2017; Paz and Guarnizo, 2019).

Extreme climatic events and increased daily climatic ranges may be more important determinants of the biological responses to climate change than temperature mean alone (Sheldon and Dillon, 2016). However, microclimates can buffer ambient temperatures and might act as refugia against such extremes (Jucker et al., 2020). This buffering could facilitate elevational shifts upwards, which could, in turn, lead to species becoming more arboreal in the cooler canopies of mountainous forests (Scheffers et al., 2014). However, highly mobile ectotherms such as flying insects often need to reach food sources hundreds of metres apart and across different forest layers, such that behavioural buffering might not be possible. Furthermore, ectotherms can have vastly different ecologies through different life stages, both in the microclimates of the forest they inhabit and in their ability to cope with thermal extremes (Klockmann et al., 2017; MacLean et al., 2016; Pincebourde and Casas, 2015; Steigenga and Fischer, 2009). Therefore, the fate of ectotherms in tropical forests will depend largely on their own thermal tolerances, as well as on the availability of local climate refugia that buffer against extreme temperatures (Nowakowski et al., 2018; Pincebourde and Casas, 2015; Scheffers et al., 2014).

Plasticity and evolutionary potential in thermal tolerance could help ectotherms cope with human-induced climate and habitat change (Hoffmann and Sgrò, 2011). Tropical species are predicted to have evolved reduced thermal plasticity compared to temperate species, due to low or absent seasonality (Sheldon et al., 2018; Tewksbury et al., 2008). A recent review found that ectotherms, in general, have low thermal tolerance plasticity, with most species having less than a 0.5°C acclimation ability in upper thermal limits (Sheldon, 2019). Detecting evolutionary change in the wild is challenging, especially in the tropics, where long-term monitoring schemes are extremely rare (Merilä and Hendry, 2014). However, in two tropical Drosophila species, moderate levels of desiccation stress have led to adaptive evolutionary responses in laboratory conditions (van Heerwaarden and Sgrò, 2014), whereas with higher and unrealistic levels of desiccation stress, flies were not able to track the changes (Hoffmann et al., 2005; Sheldon, 2019). Thus, a better knowledge of the plasticity and evolutionary potential of thermal tolerance in tropical insect species will be key to predicting their ability to cope with the warming climate, and tests with realistic levels of environmental change are required.

Accurately predicting the responses of tropical ectotherms to climate and land-use change therefore requires that we understand two complementary aspects of the system: the thermal and humidity buffering potential of tropical forests across altitudinal gradients and the thermal tolerance of the organisms inhabiting them. In this study, we: (i) measured microclimates for a full year (temperature and humidity) across the elevational range of Heliconius butterflies and assessed the accuracy of publicly available climatic predictions for the same locations, (ii) tested the heat tolerance of 10 butterfly species in the wild and (iii) reared offspring from high and low altitude populations of H. erato in common-garden conditions to test whether differences observed between wild populations were genetic or plastic.

We collected high (mean=1398 m above sea level, m.a.s.l.) and low (mean=495 m.a.s.l.) altitude populations of Heliconius butterflies with hand nets along the western and eastern slopes of the Ecuadorian Andes, at a similar latitude (Fig. S1). Every Heliconius species encountered at each site was collected, but only those species with more than five individuals at each elevation (high and/or low) were included in the analyses (15 out of 329 wild individuals removed). Detached wings were photographed dorsally and ventrally with a DSLR camera with a 100 mm macro lens in standardised conditions, and wing area was measured with an automated pipeline in the public software Fiji (following Montejo-Kovacevich et al., 2019). Butterflies with Heliconius melpomene malleti phenotypes were genotyped with a restriction digest of amplified COI genes (following Nadeau et al., 2014), to identify cryptic H. timareta spp. nov. individuals. All the images are available in the public repository Zenodo (https://zenodo.org/communities/butterfly/) and full records with data are stored in the EarthCape database (https://heliconius.ecdb.io).

Research and collecting permits were granted by the Ministerio del Ambiente, Ecuador, under the Contrato Marco MAE-DNB-CM-2017-0058.

Microclimates across elevations on both sides of the Ecuadorian Andes were characterised by recording temperature and relative humidity every hour in the understorey and mid-layers of the forest for a full year (Fig. S1), between February 2017 and February 2018. We used 40 HOBO temperature loggers (model: HOBO UA-001-08 8K; accuracy: 0.5°C, resolution: 0.1°C, 10 per area) and 16 high-accuracy humidity and temperature HOBO data loggers (model: HOBO U23-001; temperature accuracy: 0.21°C; temperature resolution: 0.02°C; relative humidity accuracy: 2.5%; relative humidity resolution: 0.05%, four loggers per area). We chose 28 forest sites that had not recently been disturbed by humans, usually inside or near nature reserves (localities in Table S1), always over 500 m away from any paved road and at least 5 m away from narrow walking trails (Movie 1). Seven of these were at high altitude (mean=1214 m.a.s.l.) and seven at low altitude (mean=444.6 m.a.s.l.) on the eastern and western slopes of the Andes (Fig. S1). Sites were at least 250 m apart from each other and in the same areas where Heliconius populations were sampled for this study. We placed one logger in the understorey (mean height=1.16 m) and one in the subcanopy (mean height=10.7 m) at each site; height was measured with a laser meter. Subcanopy loggers were as close as possible to directly above the understorey loggers, and both were hung from tree branches with fishing line and suspended mid-air (Movie 1). To prevent exposure to direct solar radiation, temperature data loggers were secured inside a white plastic bowl and humidity data loggers between two flat white plastic plates, allowing for horizontal air flow.

Heat stress resistance of wild-caught individuals was measured with a heat knockdown assay (Huey et al., 1992; Sørensen et al., 2001). Butterflies were tested less than 12 h after they were collected in the field, in the evening of the day they were collected and at approximately the same altitude. Individuals were stored in envelopes with damp cotton and fed a small amount of sugary water to standardise hydration levels before being tested. Butterflies were placed in individual glass chambers, fitted with an instant read digital thermometer (accuracy: 1.0°C, resolution: 0.1°C). Glass chambers contained 150 g of metal beads to add weight, covered with a Styrofoam platform for butterflies to stand on. Five chambers were introduced at a time into a plastic hot water bath at 51°C, and we recorded the time until the interior of the chamber, where the butterflies are placed, reached 39.0°C (‘heating up time’ hereafter). For the duration of the assay, the temperature inside the chambers was kept between 39.0°C and 41.0°C, by increasing or decreasing the hot water bath temperature as required. Tests in which the chamber temperature went below this range during the duration of the assay were removed from further analyses. We recorded the time and exact temperature at heat knockdown, defined by the loss of locomotor performance of the individual butterfly (Huey et al., 1992), i.e. when the butterfly's legs collapsed or it fell on its side. Temperature at knockout was accounted for in the models, as in 89 out of 496 assays the chamber temperature went above 41°C (Fig. S5), but assays in which temperature went above 41.9°C were removed from all analyses (n=2). Butterflies were monitored for a maximum of 60 min, and this maximum value was used for those that had not been knocked out within the time frame (n=1).

Fertilised females of H. erato were caught in the wild with hand nets in high-altitude (mean=1348 m.a.s.l.) and low-altitude (mean=543 m.a.s.l.) localities, in the vicinity of the data loggers used for microclimate analyses. Females from both altitudes were kept in separate cages of purpose-built insectaries at the Universidad Regional Amazónica Ikiam (Tena, Ecuador, 600 m.a.s.l.). Eggs were collected three times per day and individuals were reared in separate containers throughout development but stored in the same area of the insectary and randomly assorted. All offspring were individually fed the same host plant, Passiflora punctata. Development rates and pupal and adult mass were recorded for all offspring. Common-garden reared offspring were blindly tested for heat knockdown resistance in the evening 1 day after emergence following the same protocol as for wild-caught individuals. Offspring from high- and low-altitude mothers had individual IDs with no indication of their origin, thus we were able to blindly test five individuals at a time and avoid potential observer bias. Adult offspring wings were photographed and their wing areas measured with an automated pipeline in the open-access software Fiji (following Montejo-Kovacevich et al., 2019).

All analyses were run in R V2.13 (https://www.r-project.org/) and graphics were generated with the package ggplot2 (Ginestet, 2011). Packages are specified below and all R scripts can be found in the public repository Zenodo (https://doi.org/10.5281/zenodo.3634105).

Our data showed low seasonality [standard deviation of monthly averages (Bio4)=0.63], as expected for latitudes near the Equator; therefore, it was not subdivided into months. To determine the range of temperatures and humidities that butterfly populations are exposed to at different altitudes and sides of the Andes, we first estimated daily maximum, mean and minimum temperature and humidity per data logger across the year. We used linear mixed-effect models (LMMs) to determine temperature differences across forest strata (understorey and subcanopy), implemented with the lmer() function from the lme4 package (Bates et al., 2015). We included forest strata, altitude and date as fixed effects, and as random effects we placed data logger ID nested within site and area (altitude+side of the Andes). Dates were standardised to a mean of zero and unit variance to improve model convergence (Zuur et al., 2009). Significant differences between forest strata and altitudes were assessed via post hoc comparisons with Tukey tests. Mean diurnal range was estimated per day per data logger (as the daily minimum subtracted from the daily maximum) and then averaged across the year. To obtain the rate of increase in understorey temperature for every 1°C increase in subcanopy temperature, we fitted linear models with the hourly data to obtain the relationship between understorey temperature and subcanopy temperature (following González del Pliego et al., 2016).

Bio1 is more commonly known as the annual mean temperature and Bio2 as the annual mean diurnal temperature range. The latter is mathematically equivalent to calculating the temperature range for each day in a month and averaging across that month. Because our microclimate data showed low yearly seasonality (temperature s.d.=0.63) and to avoid biases in months where not all records were available (e.g. when changing logger batteries halfway through the year), we averaged across the whole year.

To test variation in thermal tolerance across species in the wild and across families in common-garden reared conditions, we first used an ANOVA approach, with species or brood as a factor explaining the variation in heat knockdown time. We then estimated within-species and within-brood trait repeatability, or intra-class correlation coefficient, with a linear mixed effect model (LMM) approach. This requires the grouping factor to be specified as a random effect, in this case species (for wild individuals) or brood (for common garden-reared offspring), with a Gaussian distribution and 1000 parametric bootstraps to quantify uncertainty, implemented with the function rptGaussian() in the rptR package (Stoffel et al., 2017). By specifying species/brood ID as a random effect, the latter approach estimates the proportion of total thermal tolerance variance accounted for by differences between species or families.

To determine the effects of altitude on thermal tolerance in wild butterflies and common-garden reared offspring, we fitted two separate LMMs, implemented with the lmer() function from the lme4 package (Bates et al., 2015). Both models initially included heat knockdown time as the response variable, and altitude at which the wild individual or mother of the brood was collected, sex, wing size (plus all two-way interactions between them), minutes that the chamber took to reach 39°C, and temperature at knockout as fixed effects. Additionally, the common-garden reared offspring model included development time (days from larva hatching to pupation) and brood egg number (to control for time the brood mother spent in the insectary) as fixed effects. The random effects for the wild individual model and for the common-garden reared offspring model were species identity and brood mother identity, respectively. All fixed effects were standardised to a mean of zero and unit variance to improve model convergence (Zuur et al., 2009).

We implemented automatic model selection with the step() function of the lmerTest package (Kuznetsova et al., 2017), which performs backwards selection with likelihood ratio tests and a significance level of α=0.05, first on the random effects and then on the fixed effects to determine the simplest best-fitting model. Model residuals were checked for homoscedasticity and normality. To obtain P-values for the fixed effects, we used the anova() function from the lmerTest package, which uses Satterthwaite approximation (Kuznetsova et al., 2017). We estimated with the coefficient of determination (R²) the proportion of variance explained by the fixed factors alone (marginal R², R_LMM,m²) and by both the fixed effects and the random factors (conditional R², R_LMM,c²), implemented with the MuMIn library (https://CRAN.R-project.org/package=MuMIn; Nakagawa et al., 2017).

Overall, the patterns of differentiation between forest layers and altitudes are very similar across sides of the Andes in our study (Fig. 1A versus B), and there was low seasonality across months (Fig. 1). Our measured temperature seasonality per logger, averaged across areas, was 0.61 and 0.81 for the western highlands and lowlands, respectively, and 0.49 and 0.60 for the eastern highlands and lowlands, respectively. The lowland sites were, on average, 4.1°C warmer than the highland sites, which were over 750 m apart in elevation (Tables S1). Annual mean temperatures interpolated from WorldClim2 were always closer to subcanopy strata annual means (Fig. S2A). In contrast, the understorey annual mean temperatures were on average 0.5°C lower than WorldClim2 annual mean temperatures (dotted lines in Fig. 1Ai,ii,Bi,ii, Table S1). The minimum temperatures were consistently higher in our microclimate data compared with the interpolated monthly minima estimated in WorldClim2, especially for the highlands (lower dotted line in Fig. 1Aii,Bii).

Forest canopies thermally buffered the understorey, but more so during the day and in the lowlands. During the day, the understorey thermally buffered the canopy temperature maxima by 1.19–1.62°C in the highlands, and by 1.98–2.11°C in the lowlands (Table 1). At night, understories buffered the temperature minima of the canopies by 0.07–0.12°C in the highlands and by 0.25–0.23°C in the lowlands (Table 1). Thus, the forest buffered high and low temperatures in the understorey throughout the day and night, respectively (Fig. 2A). Temperature differences between day and night are greater in the lowlands, where days are warmer (Fig. 2A), but less so in the understorey of all areas. On average, the difference between subcanopy and understorey diurnal thermal range in the highlands was 1.34°C, whereas in the lowlands this difference was 2.09°C (Fig. 2B). However, WorldClim2 interpolations for diurnal thermal ranges were 3.5°C higher than our records in the highlands, resulting in the highlands being predicted to have higher thermal ranges than the lowlands (stars, Fig. 2B). This was the opposite elevational trend to that observed in our data, where thermal ranges were lower in the highlands.

The temperature buffering of the lowlands was higher than the highlands, so that for every 1°C increase in subcanopy temperature the understorey increased by 0.68°C in the lowlands, in contrast to 0.73°C in the highlands. The lowland canopies exceeded 39°C on 31 days throughout the year, whereas the highlands never did (Fig. 1). The monthly maximum temperature interpolated by WorldClim2 was close to our measured subcanopy temperature maxima in the lowlands, but 2.01°C higher, on average, in the highlands (Fig. 2C). WorldClim2 monthly minima were also overestimated at both elevations, predicting it to be 2.6°C cooler at night in the highlands and 2.3°C cooler in the lowlands than our measured understorey microclimates (Fig. S2B).

In the highlands, the understorey daily relative humidity minimum was on average 3.7 percentage points higher than in the subcanopy, whereas in the lowlands the difference was 11.8 percentage points (Fig. S3). The drying power of the air (VPD) varied drastically between layers of the forest, but more so in the lowlands (Fig. 3). In the highlands, maximum daily VPD from the same site averaged 1.1 hPa higher in the subcanopy compared with the understorey, whereas in the lowlands the difference was 4.8 hPa. The threshold for tree transpiration is thought to be at 12 hPa in the tropical montane areas, above which transpiration, and thus growth, is impeded by the drying power of the air (Motzer et al., 2005). This threshold was exceeded 883 times in the lowland subcanopy across our data loggers, 102 times in the lowland understorey, 12 times in the highland subcanopy and never in the highland understorey (Fig. S3).

Heat tolerance varied across species (ANOVA: F_9,268=8.75, P<0.0001; Fig. 4A), with 44% of this variation explained by species identity (repeatability=0.44, s.e.=0.13, P<0.0001). In the highlands, 33% of the individuals tested were knocked out before the chambers reached a temperature of 39°C, in contrast to 7% of the individuals tested in the lowlands (n_high=60/183, n_low=6/95; Fig. 4B). Mean (±s.e.m.) heat tolerance across species was on average 5.4±0.57 min for highland individuals and 15.9±1.43 min for lowland individuals (red dashed line in Fig. 4A). Altitude and time until the chamber reached 39°C were significant predictors of knockout time in wild individuals (Table S1), with the fixed effects alone explaining 29% of the variation in thermal tolerance (R_LMM,m²=0.29) and 39% when considered together with species identity as a random effect (R_LMM,c²=0.39). The high-altitude populations of the two most evenly sampled species across altitudes, H. erato and H. timareta, were less thermally tolerant than their lowland conspecifics (t-test, H. erato: t₇₇=−5.3, P<0.0001, H. timareta: t₁₄=−2.3, P<0.05; Fig. S6). In this part of the eastern Andes of Ecuador, the lowland H. melpomene have been largely replaced by a cryptic subspecies of H. timareta (Nadeau et al., 2014), which explains the low numbers of lowland individuals of H. melpomene, a species with a very wide range across the Neotropics.

Common-garden reared offspring of H. erato lativitta varied in heat tolerance across families (ANOVA, F_14,262=5.15, P<0.0001), and 25% of this variation was explained by brood identity (repeatability=0.25, s.e.=0.10, P<0.0001). In the wild, low-altitude H. erato lativitta were, on average, able to withstand high temperatures for 10 min longer compared with high-altitude populations (t-test: t₇₇=−5.3, P<0.0001; Fig. 4, left). In contrast, when reared in common-garden conditions, individuals from lowland broods were able to withstand heat for only 1.4 min longer than offspring from highland broods (Fig. 5, right). As a consequence, parental altitude only had a marginally significant effect on offspring thermal tolerance (knockdown time), whereas experimental variables, such as time until the chamber reached 39°C and the temperature at knockout, were significant predictors of knockout time in the offspring (Table 2, ‘common-garden reared model’). Fixed effects alone explained 39% of the variation in thermal tolerance (R_LMM,m²=0.39) and 48% when together with brood identity as the random effect (R_LMM,c²=0.48), indicating trait heritability (Table 2). The variance in thermal tolerance in wild populations was higher than in common-garden reared offspring (Fig. 5), probably owing to higher variation of developmental temperatures, fitness and age in the wild. The number of eggs a mother had previously laid had a positive and significant effect on adult thermal tolerance, and interacted with parental altitude, likely owing to high-altitude mothers living longer in the insectary.

We found that tropical forests have great climatic buffering potential, especially at lower elevations, and that this was similar across two independent elevational clines on both sides of the Ecuadorian Andes (Figs 1 and 2). Interpolated climatic variables for these same areas did not capture our observed microclimates, especially at high elevations (∼1100 m.a.s.l.), where weather stations are very sparse (Fick and Hijmans, 2017). Furthermore, we found evidence for differences in thermal tolerance in the wild across 10 butterfly species, regardless of whether they had altitudinally narrow or widespread distributions (Fig. 4). However, these differences were greatly reduced when a widespread species was reared in common-garden conditions, indicating likely non-genetic, environmental effects on temperature tolerance.

Many macroecological predictions of species distributions have emerged from the notion that tropical latitudes lack strong seasonality (Janzen, 1967; Sheldon et al., 2018; Shelomi, 2012). However, while generally true, these predictions ignore the climate complexity of tropical forests. Our results show that lowland and montane tropical forests buffer ambient temperature and humidity. The understory daily maximum temperature was, on average, 1.4°C and 2.1°C cooler than the subcanopy, at high and low elevations, respectively (Table 1). The temperature offset between forest canopy and understory becomes larger at more extreme temperatures (De Frenne et al., 2019). Thus, lowland forests buffered temperatures to a greater extent than highland forests, which is of particular importance for ectotherms in lowland environments, which are routinely exposed to extreme temperatures (Deutsch et al., 2008). It is important to note that in the present study, the highland areas were located at mid-elevations and our lowlands were in the foothills of the Andes, which have been less studied because the differences are often assumed to be marginal. Similarly, our subcanopy data loggers were positioned 10 m from the forest floor, thus we would expect the temperature and humidity offsets to be even stronger when compared with above-canopy or open-habitat temperatures.

Janzen's (1967) hypothesis predicted that the reduced climatic variability in the tropics would result in ectotherms having narrower thermal tolerances, and, in turn, reduce dispersal across elevations (Sheldon et al., 2018). In the present study, the mean temperature difference between the subcanopy and the understorey, only 10 m apart, was 0.25°C in the highlands and 0.44°C in the lowlands (Table 1). This is more than the temperature change across these elevational clines, where for every 10 m in elevation there was a 0.05°C decrease in temperature. Yet, in the wild, we found that low-elevation populations were much more tolerant to heat than highland populations. Although in this study we did not measure cold tolerance, we found that the difference in minimum temperatures across elevations and forest strata is much smaller than that of maximum temperatures (Fig. 1Av,Bv). Thus, we can hypothesise that, given a linear change in cold tolerance, the disproportionately greater heat tolerance of lowland Heliconius would result in them having broader thermal breadths than high-elevation populations. The microclimatic variability that tropical ectotherms are exposed to within their habitats might offset the lack of seasonality across the year, making some species more able to cope with warming than others. Thus, protecting tropical forests' climatic buffering potential across elevations is essential to enable potential upland shifting in the face of climate and land-use change.

We found a clear disparity between field-collected microclimate and interpolated macroclimate temperature (Fig. 2). For instance, WorldClim2 estimated the maximum daily temperatures to be 2.01°C higher in the highlands, and overestimated annual mean temperature in the understorey by 0.5°C. Strikingly, these values are similar in magnitude to the projected warming for the next century in Andean mid-elevations (Beaumont et al., 2011; Urrutia and Vuille, 2009). In large part, these disparities can be attributed to the fact that coarse-gridded temperature surfaces, such as WorldClim2, are interpolated from weather stations that are located in open habitats (De Frenne et al., 2019; De Frenne and Verheyen, 2016). Several recent studies have reported similarly striking differences (Jucker et al., 2018; Lembrechts et al., 2019; Potter et al., 2013; Storlie et al., 2014). Near-surface temperatures can only be accurately measured with in situ loggers or with emerging remote sensing technologies, such as airborne laser scanning (Jucker et al., 2018). Furthermore, WorldClim2 interpolations at high altitudes tend to be less accurate, especially in the tropics, where weather station data are very sparse (Fick and Hijmans, 2017). This raises the question of how useful coarse macroclimatic grids are for assessing thermal tolerances of organisms that are affected by fine-scale microclimates (De Frenne et al., 2019; Lenoir et al., 2017; Navas et al., 2017; Nowakowski et al., 2018). In addition, very few studies in the tropics have accounted for humidity and VPD variability at the microclimate level (Bujan et al., 2016; Friedman et al., 2019; García-Robledo et al., 2016), as the loggers required to do so can be four to five times more costly than temperature loggers. Nevertheless, inclusion of VPD in species distribution models has been shown to significantly improve their accuracy (de la Vega and Schilman, 2017). The differences in VPD between subcanopy and understorey observed here were much more pronounced in the lowlands, highlighting the importance of protecting forest complexity in these areas, which are under constant threat of land-use change.

Our extremely limited knowledge of thermal tolerance and plastic potential of tropical ectotherms in the wild further hinders our ability to predict ectotherm responses to climate change. As expected, the ability of wild butterfly species to cope with extreme, but natural, levels of heat (∼40°C) was much lower for those inhabiting high altitudes. The heat tolerance differences between highland and lowland populations of widespread species were much lower than between highland and lowland specialist species (Fig. 4, Fig. S6). This suggests that although plasticity plays a role in widespread species, as shown in our common-garden rearing with H. erato, elevationally restricted species are likely to have fixed genetic characteristics that make them better suited to their thermal environment, and thus show more differences when compared across elevations. Behavioural shifts might alleviate the impact of climate extremes in tropical forests, but the capacity of a species to shift sufficiently will be constrained by life history, energy budgets and thermal tolerance (Sheldon and Tewksbury, 2014). In these butterflies, altitude has been shown to pose strong selective pressures, and some are constrained in their body size by contrasting reproductive strategies (gregarious versus solitary), which could, in turn, restrict adaptive plastic responses to environmental change (Montejo-Kovacevich et al., 2019). Nevertheless, the observed thermal buffering of forests would undoubtedly benefit and be exploited by these butterflies during extreme temperature events, such as the 31 days where temperatures in the lowlands went above 39°C.

When reared under a common developmental temperature, the differences in adult thermal tolerance observed in the wild largely disappeared in the widespread species, H. erato lativitta, indicating plasticity in individual thermal tolerance (Fig. 4). In this study, we cannot distinguish developmental plasticity from acclimation. The former would imply that larval rearing temperature would irreversibly determine adult heat knockdown resistance, whereas the latter implies a reversible response to temperature (Llewelyn et al., 2018). Interestingly, we did find evidence of at least some genetic component to heat knockdown resistance, as 25% of its variation among common-reared offspring was explained by family identity. Thus, it is likely that a combination of genetic and environmental effects determine adult thermal tolerance in Heliconius butterflies. Heat knockdown resistance has been associated with the expression of a heat shock protein (Hsp70) and its plasticity has been shown to vary across altitudinally and latitudinally structured populations of Drosophila buzzati and European butterflies (Karl et al., 2009; Luo et al., 2014; Sørensen et al., 2001), suggesting a mechanism for both genetic and environmental control of heat knockdown resistance. Highly heritable thermal performance traits have been shown to allow rapid adaptation to changing climates in lizards (Logan et al., 2014), but developmentally plastic traits with varying and unknown levels of heritability have less predictable roles (Merilä and Hendry, 2014).

The fact that 45% of the individuals from three high-altitude Heliconius species were knocked out at temperatures between 35 and 39°C is remarkable (Fig. 3B), given that in the highlands, there were 50 days in a year where temperatures went above 30°C. These high-altitude species would have suffered high mortality if reared in our lowland insectaries, as temperatures often rise above 35°C. In contrast, only 27% of sympatric high-altitude individuals of widespread species were knocked out at temperatures between 35 and 39°C. Canopy-dwelling tropical ants are known to behaviourally circumvent high-temperature areas encountered while foraging, avoiding temperatures 5–10°C below their thermal tolerance limits (Logan et al., 2019; Spicer et al., 2017). Thus, these butterflies and other high-altitude species must adapt their flying times and/or behaviours during hot periods, which could have cascading effects on fitness. Heliconius follow the same flower foraging trap-lines every day, as they require large amounts of nectar and, uniquely for lepidopterans, pollen, to thrive throughout their long adult life-spans (Jiggins, 2016; Mallet, 1986). Hot patches of the forest close to their thermal limits could severely hinder their ability to follow a foraging path, disrupting a basic biological function. Thus, the climatic refugia that the understorey provides could be crucial to cope with the ongoing increase in extreme temperature events in high elevation habitats (Ruiz et al., 2012; Scheffers et al., 2014).

Relative humidity and VPD, which have received little attention in the literature, were also greatly buffered by forest canopies (Fig. 3). High VPD has direct impacts on seedling growth and survival (Jucker et al., 2018; Motzer et al., 2005), as well as on ectotherm activity levels through increased desiccation risk (Bujan et al., 2016; Friedman et al., 2019). In canopy tropical ants, desiccation resistance was negatively correlated with thermal tolerance, suggesting an energetic trade-off between the two traits, but both were generally higher than in understorey ants (Bujan et al., 2016). Thus, forest temperature and humidity buffering may benefit canopy-dwelling butterflies disproportionately. Furthermore, understorey specialist species are thought to be closer to their thermal limits owing to the stable climatic conditions they inhabit, and thus these may not have suitable refugia or the thermal capacity to cope with warming (Huey et al., 2009; Scheffers et al., 2017; Sheldon, 2019). Heliconius erato is considered an open-habitat dwelling species, whereas others not studied here, such as H. sapho and H. cydno prefer closed-canopy forest (Estrada and Jiggins, 2002). Future research could focus on testing desiccation risk across Heliconius species with different flying habits and at different elevations. Increased detrimental effects on growth and survival under climate change could also hold true for less mobile, understorey-dwelling life stages, such as butterfly larvae.

Habitat degradation and land-use change in the Andean foothills are pressing concerns for tropical conservation. Ectotherms escaping unsuitable climates in the lowlands may struggle to shift their ranges upwards owing to habitat fragmentation (Chen et al., 2009; Scheffers et al., 2016). For example, most remnants of pristine forest in Borneo are in montane areas, as the flat lowlands have been converted to oil palm plantations. Current protected areas may fail to act as stepping stones as they are too isolated and distant from upland climate refugia (Scriven et al., 2015). A similar scenario occurs in western Ecuador, one of our study areas, where a large portion of the low and mid-elevations have been converted to agricultural lands in the past three decades (Wasserstrom and Southgate, 2013). Any habitat change that affects forest heterogeneity could reduce its large temperature buffering potential (Blonder et al., 2018; Jucker et al., 2020), and butterfly diversity as a whole (Montejo-Kovacevich et al., 2018). Nevertheless, microclimates have been shown to recover decades after low impact land uses (González del Pliego et al., 2016; Mollinari et al., 2019; Senior et al., 2018), allowing for recolonization of biodiversity (Hethcoat et al., 2019). This highlights the need to protect degraded secondary forest, as these are now more abundant than primary forests in most of the tropics (Edwards et al., 2011; Senior et al., 2017).

Tropical ectotherms find themselves in highly heterogeneous, threatened habitats, which have greater climate buffering potential than previously thought (De Frenne et al., 2019; De Frenne and Verheyen, 2016). However, the low seasonality of tropical environments together with the steep environmental gradients of montane habitats, as Janzen's hypothesis predicts, makes tropical ectotherms particularly vulnerable to climate and habitat change (Deutsch et al., 2008; Huey et al., 2009; Janzen, 1967; Polato et al., 2018). The clear mismatch between the microclimates we measured and the widely used interpolated global datasets highlights the importance of field-based climate measurements. Furthermore, the striking difference in heat tolerance and evidence for plasticity across a moderate elevational cline demonstrates the importance of temperature for the persistence of tropical ectotherms. Our results suggest that the inclusion of microclimate buffering within models and experimental testing of thermal tolerances is crucial for incorporating realistic temperatures experienced by small organisms (Paz and Guarnizo, 2019). More research into the evolvability and plasticity of heat tolerance is needed to accurately assess the vulnerability of tropical ectotherms in the face of anthropogenic change.

We would like to thank all field assistants that have contributed to this study, particularly to Ismael Aldas, Jennifer Smith, and Melanie Brien. We are extremely grateful to the Mashpi Reserve, Narupa Reserve (Jocotoco Foundation), Jatun Satcha Reserve, and Universidad Regional Amaóznica Ikiam for their support.