Colour variation across climatic gradients is a common ecogeographical pattern; yet there is long-standing contention over underlying causes, particularly selection for thermal benefits. We tested the evolutionary association between climate gradients and reflectance of near-infrared (NIR) wavelengths, which influence heat gain but are not visible to animals. We measured ultraviolet (UVA), visible (Vis) and NIR reflectance from calibrated images of 372 butterfly specimens from 60 populations (49 species, five families) spanning the Australian continent. Consistent with selection for thermal benefits, the association between climate and reflectance was stronger for NIR than UVA–Vis wavelengths. Furthermore, climate predicted reflectance of the thorax and basal wing, which are critical to thermoregulation; but it did not predict reflectance of the entire wing, which has a variable role in thermoregulation depending on basking behaviour. These results provide evidence that selection for thermal benefits has shaped the reflectance properties of butterflies.

1. Introduction

Variation in animal coloration forms the basis of some of the earliest described and best-known ecogeographical patterns. Such patterns include Bogert's rule (darker individuals in cooler environments [1,2]) and Gloger's rule (darker individuals in warmer and more humid (tropical) than cooler, dryer (temperate) environments [3,4]). A commonly invoked explanation for these ecogeographical patterns is selection for thermal benefits; but this explanation remains contentious for several reasons [5,6]. First, colour is influenced by selection for multiple competing functions (e.g. camouflage, signalling, protection from ultraviolet and bacterial damage [7]) and different selection pressures can produce similar ecogeographical patterns [8]. Second, the relationship between colour and heat gain is complex and depends strongly on environmental conditions, which vary considerably among microhabitats and throughout the day and seasons [9]. Thermal benefits and costs of coloration can also be enhanced or compensated by behaviour, such as postural changes or seeking shade [5]. Third, most studies have only considered reflectance properties over human-visible wavelengths (Vis, 400–700 nm), sometimes in conjunction with the ultraviolet range (UVA, 300–400 nm) [5,6]. Few have considered reflectance properties beyond 700 nm. This is important because wavelengths of direct sunlight greater than 700 nm (the near-infrared, NIR) comprise about approximately 50% of total incident radiation (electronic supplementary material, figure S1) [5]. Consequently, visible colour can be a poor predictor of an organism's absorption of solar radiation [9,10].

A powerful way to test whether selection for thermal benefits has shaped animal reflectance properties is to assess the evolutionary association between climate and reflectance in unseen NIR wavelengths. We expect NIR variation to be primarily shaped by selection for thermal functions because NIR wavelengths strongly affect heat gain but cannot be seen by most animals [5]. By contrast, we expect the relationship between visible reflectance (colour) and climate to be weak or variable because visible reflectance, particularly of the wings, is more likely to be constrained by selection for optical functions such as camouflage, sexual signalling and warning coloration. Comparison of drivers of visible and NIR reflectance variation also provides a means to test whether animals have ‘solved’ the problem of competing selection for optical and thermal functions by modulating visible and NIR reflectance differently. For example, low visible reflectance for camouflage (e.g. brown coloration) may be offset by high relative NIR reflectance to mitigate the risk of overheating. Recently, we showed that birds in hot, dry climates have higher relative NIR reflectance than species occupying more thermally benign environments, irrespective of visible colour [11]. This pattern was most evident for smaller species, which gain greater benefits from higher reflectivity through reduced evaporative water loss. However, the thermal significance of NIR reflectivity is expected to be even greater for ectotherms, which lack the buffering effect of a complex insulating layer (fur or feathers).

Butterflies are excellent models for testing the thermal benefits of reflectance and resulting ecogeographical patterns. Butterflies are primarily ectothermic and regulate their temperature by basking in sunlight. This includes dorsal basking with wings open to directly expose the body and dorsal surfaces to the sun, and ventral basking whereby the wings are folded to expose the ventral face to the sun [12]. The wings are important for regulating thoracic temperature because butterflies can close their wings to shade the thorax or adjust the angle of the wings to reflect solar radiation onto the thorax [13]. However, the thorax and basal region of the wing, the first few millimetres adjacent to the thorax, are the most important body regions for thermoregulation because they directly absorb and conduct heat to the flight muscles [14]. There is little fluid and vascular extension through the wings beyond the basal region and the rate of haemolymph flow is too low to transfer significant quantities of heat to the thorax via circulation [14–16]. Given the roles of the different body regions in thermoregulation, we expect a more consistent relationship between climate and the reflectance properties of the thorax and basal wing than the entire wing.

Here, we use a phylogenetic comparative approach and digital imaging to assess the relationship between climate and reflectance properties from 320 to 1100 nm. This range represents the limits of digital camera sensors but captures approximately 83% of energy in incident sunlight. We used calibrated images to quantify UVA–Vis (320–700 nm) and NIR (700–1100 nm) reflectance on the thorax, basal wing region and entire wing of 372 specimens of butterflies (60 colour races, 49 species, five families) spanning the range of climates across the Australian continent. Specifically, we tested the predictions that (i) the association between climate and NIR reflectance is stronger than the corresponding association with UVA–Vis reflectance, which may have functions other than thermoregulation; and (ii) that the association between climate and reflectance properties of the thorax and basal wing is stronger than the corresponding association with the entire wing, which has a variable role in thermoregulation. We find support for these predictions, providing evidence that selection for thermal benefits has shaped the reflectance properties of butterflies.

2. Material and methods

(a) Specimen sampling

Butterfly specimens were sampled from the entomology collection at Melbourne Museum, Australia. We photographed specimens from the five butterfly families with non-migratory species in Australia (Nymphalidae, Lycaenidae, Hesperiidae, Pieridae and Papilionidae) to represent a broad diversity in phylogeny, colour, size, sexual dimorphism, distribution and climatic diversity. We sampled distinct colour populations for species which comprised more than one geographical colour form. For species that comprised a single widespread colour form in the collection, we sampled from the geographical location that had the highest number of best-preserved specimens. We photographed 12 colour populations within each family, totalling 60 populations from 49 species. We photographed three to five individuals for the 43 populations where we could not identify the sex of the specimens from obvious external features, and for the remaining 17 populations with distinguishable sexes, we photographed three to five individuals per sex, totalling 372 specimens (electronic supplementary material, table S1). Our specimen selection was restricted to colour populations that had a minimum of three well-preserved specimens with wings opened and angled relatively flat, and minimal structural damage to their body and wings. Neither UVA–Vis nor NIR reflectance were correlated with specimen age (electronic supplementary material).

(b) Imaging set-up

We used a Nikon D7200 digital SLR camera with a full spectrum fused-silica conversion (Camera Clinic, Melbourne) and a high performance 60 mm UV–Vis-IR Apo Macro 1 : 4 lens with quartz optics (CoastalOpt, JenOptik, USA). The lens was fitted with different optical filter combinations (Edmund Optics, Singapore; electronic supplementary material figure S2) to enable photography of the specimens within three spectral bands: UVA (320–400 nm), Vis (400–700 nm) and NIR (700–1100 nm). A modified flash unit (Nikon SB-140 UVIR clone; BeyondVisible, USA) was used as a light source for UVA and Vis photos (2.5 mm open-slit filter for Vis images; spectral power distribution provided in electronic supplementary material, figure S3), and a 150 W halogen light (Arlec HL108) with a 78 mm linear tungsten halogen tube (HT150) was used for NIR photos. We also used a custom-made, parabola-shaped aluminium reflector mounted on the opposite side of the camera (electronic supplementary material, figure S4) to evenly illuminate the subjects.

Two different camera heights, and therefore magnification ratios (1 : 4 and 1 : 5), were used to photograph small and large butterflies. Photography conditions and settings were determined based on extensive testing of temporal stability and spatial consistency in illumination, fully detailed in electronic supplementary material, figures S5–S7 and table S2). We used photography rather than spectrophotometry to enable us to measure reflectance over the entire body regions of interest.

(c) Butterfly imaging

All 372 butterflies were individually photographed in a dark room three times to generate a UVA, Vis and NIR image (figure 1). Each image included a Spectralon reflectance standard (LabSphere, NH, USA) reflecting 42% of all incident radiation within a 320–1100 nm spectral interval, which also served as a scale bar. Each specimen was pin-mounted in a piece of grey foam with the thorax aligned at the height of the 40% grey standard (lens focus-point). The resulting 1116 images were batch encoded into 300 dpi, 8-bit Adobe 1998 RGB colour space [17] TIFF files using Adobe Camera Raw (9.7.0) [18] for subsequent processing. Figure 1. Left to right: UVA, visible and NIR calibrated photographs of (a) Eurema hecabe phoebus (unknown sex), (b) Arhopala wildei wildei (female), (c) Graphium sarpedon choredon (unknown sex) and (d) Ogyris genoveva genua (male). Lighter colours (higher intensity) indicate higher reflectance in calibrated images. Note differences in intensity and contrast between colour pattern elements in the UVA, visible and NIR photos.

To quantify reflectance, we extracted red (R), green (G) and blue (B) pixel intensity values from the 40% grey standard, thorax, basal wing region and entire wing area for all photos using the built-in ‘RGB Measure’ plugin in ImageJ v. 1.50i [19]. Pixel intensity values were recovered from image selections delimited with the freehand selection tool (electronic supplementary material, figure S8). Images were scaled with the 3 cm diameter of the grey standard in each image to measure the physical dimensions (wing area) using the scale tool available in ImageJ. The stacked images could then be swapped while maintaining the highlighted section, allowing the RGB pixel values to be extracted from the exact same location in each image per specimen.

We used the extracted RGB pixel intensity values to estimate mean reflectance for each of the UVA, Vis and NIR wavebands. For each waveband, we used data from one or more channels depending on the sensitivity of the camera's three sensors (B channel for UVA; R, G and B channels for Vis; R and B channels for NIR; electronic supplementary material). We calibrated RGB values based on the camera's responses to six grey standards with known reflectance (SphereOptics; Herrsching, Germany; 5%, 10%, 25%, 50%, 80% and 99% standards; electronic supplementary material, figure S9) [20,21]. Specifically, we derived a bi-exponential linearization function of the form y = aebx + cedx, where y is the linearized pixel value and a, b, c and d are empirically derived constants specific for each colour channel of our camera [22]. We derived a separate linearization equation for each selected channel (R, G and B) at each camera magnification setting (electronic supplementary material, table S3). We equalized the linearized RGB values relative to the RGB pixel values of the grey reflectance standard in each image [20,21]. We then averaged the selected channels for each waveband to get a single average reflectance value for each waveband (UVA, Vis and NIR). We also calculated the reflectance (%) for the UVA–Vis waveband combined (320–700 nm), and total reflectance (UVA–NIR: 320–1100 nm). Specifically, the UVA waveband is 320–400 nm = 80 nm and the Vis is 400–700 nm = 300 nm, so the mean reflectance for the UVA–Vis combined = [(a80 + b300)/380], where a and b are the mean reflectance from UVA and Vis photos, respectively.

Assuming direct normal standard irradiance for dry air (ASTM G173–03) derived from SMARTS (2.9.2) [23], the wavelength range 320–1100 nm represented by our images incorporate a total of 83.27% of incident solar radiation, comprising UVA (320–400 nm; 6.7%), Vis (400–700 nm; 42.2%) and NIR (700–1100 nm; 34.4%). Thus, less than 17% of energy in sunlight falls at NIR wavelengths beyond 1100 nm and is not captured by our images. Reflectance properties of biological materials are much less variable beyond about 1400 nm [10,11,24], so our images capture most of the biologically meaningful variation in NIR reflectance.

(d) Phylogeny

We used the eight loci (COI, EF-1a, Wingless, RpS5, MDH, GAPDH, CAD and IDHT) dataset of Heikkilä et al. [25] as the backbone to our phylogenetic analyses. Where possible, we extended this dataset with corresponding published sequences for the focal species and/or congenerics [26–31] (electronic supplementary material, table S4). The resulting final backbone dataset contained 126 species, and to this, we added an additional 35 taxa for which only COI sequences were available. Each gene was aligned separately using the global-pair (G-INS-i) algorithm in MAFFT (v. 7.310) [32]. The best model of sequence evolution for each gene was identified using ModelFinder as implemented in IQTREE [33]. There were 11 cases where sequences were not available for subspecies or colour populations examined in this study. In these cases, we simulated COI sequences with the program Seq-Gen [34]. Specifically, available conspecific sequences were used as the initial state, with substitutions simulated given a model of sequence evolution (best model chosen and parameters estimated using IQTREE) and an expected final level of divergence from the initial state of 1.5%. A time-calibrated, phylogenetic estimate was then obtained using BEAST2 [35,36] with three fossil calibrations as per Heikkilä et al. [25], a relaxed lognormal model and the Yule tree prior. Two independent runs of 50 million generations (logged every 10 000) were executed and subsequently assessed for convergence. The maximum clade credibility (MCC) tree for the full dataset of 161 taxa is presented in electronic supplementary material, figure S10. A final subset of 4500 post-burnin trees (subsampled using logcombiner) were then pruned of all non-focal taxa using the R package Phylotools [37]. These trees were used in subsequent phylogenetic comparative analyses (figure 2). Figure 2. Phylogenetic relationships (MCC tree) across species/populations and families for specimens sampled (49 species, 60 populations) and average values of UVA–Vis and NIR reflectance per species in specimens sampled, for three different body regions. (Online version in colour.)

(e) Climate data

Climate data were compiled for each specimen locality using the 5 km (0.05°) grid resolution Australian Bureau of Meteorology Australian Water Availability Project (AWAP) [38] database, interpolated from weather stations. We used data from the 1 January 1990, the first year in which solar radiation data became available, to the 31 December 2016. We extracted daily data nearest to each specimen's locality for the following climate variables: rainfall (millimetre), maximum and minimum air temperature (°C), vapour pressure (hPa) and solar radiation measured as daily radiant exposure (MJ m−² d−1). Cloud coverage (%) was also calculated from the ratio of observed integrated daily solar radiation to that of a cloudless day, estimated from the NicheMapR microclimate model [39]. Using the daily values from 1990 to 2016, we calculated a single monthly mean value of each variable corresponding to each specimen over the months of activity given for the relevant population or species in The complete field guide to butterflies of Australia [40]. For example, a colour population had five monthly average values per variable if it had an activity window of five months at the specimen's collection locality. We then only chose the lowest monthly average value (average minimum value, termed ‘low’) and the highest monthly average value (average maximum value, termed ‘high’) for each variable to represent the extreme macro-climate ranges specific to each colour population and specimen's location and active months.

We performed a principal components analysis on a scaled data matrix of eight climate variables with pairwise correlation coefficients of r < 0.8: average monthly rainfall in the driest month and wettest month of activity, average monthly day-time temperature in the coldest month and hottest month of activity, average monthly solar radiation in the least radiative month and most radiative month of activity, and average monthly cloud coverage in the sunniest month and cloudiest month of activity.

Principal components (PCs) 1 and 2 had eigenvalues greater than 1 and accounted for 70.9% of the variation (47.8% and 23.1%, respectively, electronic supplementary material, table S5). PC1 loaded most strongly and positively against average monthly rainfall in the wettest month of activity, and average monthly day-time temperature in the coldest month and driest month of activity. Higher values of PC1 represent warmer, tropical climates with PC1 increasing towards lower latitudes (figure 3). PC2 loaded positively against average monthly rainfall in the driest month and average monthly cloud coverage in the cloudiest month, and negatively against average monthly solar radiation in the most radiative month. High values of PC2 represent wetter, cloudy climates, with the highest values of PC2 found in montane areas from the Wet Tropics and Tasmania (figure 3). Figure 3. Geographical plots of climate PC1 and PC2 values for each specimen at their collection coordinates (n = 372). The size of the points is proportional to the magnitude of the PC value. Map colours show the average temperature across Australia, with redder values representing hotter environments (taken from www.ala.org.au).

(f) Phylogenetic comparative analysis

To assess variation in NIR reflectance across species, accounting for variation in UVA–Vis reflectance, we calculated relative NIR reflectance as the average species residuals for a given level of UVA–Vis reflectance from a phylogenetically corrected linear regression between log 10 UVA–Vis and NIR reflectance. We used log 10 UVA–Vis reflectance to produce a linear relationship with NIR reflectance (electronic supplementary material, figure S11) because the relationship between UVA–Vis and NIR reflectance is nonlinear [11,24].

Next, we tested whether dorsal reflectance on the thorax, basal wing region and entire wing was predicted by climate and wing size using phylogenetic generalized least-squares regression (PGLS) implemented in the package ‘caper’ [41] in R v. 3.4.1 [42]. We accounted for the effects of phylogeny on trait similarities between the sampled colour populations by running each analysis on 4500 trees from the posterior distribution obtained in the Bayesian analysis. For each body region (thorax, basal wing and entire wing), we ran models with UVA–Vis reflectance and NIR reflectance as response variables. The predictors were PC1, PC2, log 10 wing area, and the interaction between PC1/ PC2 and log 10 wing area. We also included log 10 UVA–Vis reflectance as a covariate in linear models predicting NIR reflectance so we could gauge the relationship between NIR reflectance and climate, accounting for variation in UVA–Vis reflectance. We used the command dredge in the package ‘MuMIn’ to compare all possible models and select the best model according to the AICc value [43]. We plotted the distribution of residuals to assess model fit and convergence. For each model we report parameters for each predictor in the best model in the majority of trees, we also report the percentage of trees for which model was the best one and the average evidence ratio (ER) between the best model and the null model (when UVA–Vis reflectance was the response) or the model with only log 10 UVA–Vis as predictor (when NIR was the response variable). The ER provides a measure of how much more likely the best model is than the comparison model [44]. To graph our results, we calculated the MCC tree across the 4500 trees from the posterior distribution, and ran reduced models using this tree. Trend lines in graphs were predicted using these reduced models. We built the graphs presented using the packages ‘Rmisc’ [45], ‘ggtree’ [46] and ‘ggplot2’ (v. 1.5) [47].

3. Results

(a) Variation in reflectance

The UVA–Vis and NIR reflectance were highly correlated (correlation between log 10 UVA–Vis and NIR reflectance: p < 0.001, R2 = 0.77; electronic supplementary material, figure S11). Nevertheless, there was substantial variation in NIR reflectance across species (figure 2), even after accounting for visible reflectance. Relative NIR (calculated as the average species residuals for a given level of UVA–Vis reflectance from a phylogenetically corrected linear regression between log 10 UVA–Vis and NIR reflectance) ranged from very low (−12.603) in Cressida cressida to very high (11.06) in Paralucia pyrodiscus, and indicated up to around 20% variation in NIR for a given level of UVA–Vis reflectance. Variation in NIR within species was very low, the mean SD within species was 2.73% and 90% of the species had SD below 5% reflectance. For the subset of species for which the sexes could be easily distinguished, there were significant differences in reflectance between sexes, with females having higher reflectance (UVA–Vis and NIR) of the thorax and basal wing, and males having higher reflectance (UVA–Vis only) of the entire wing (electronic supplementary material, figure S12). Phylogenetic signal was very high for UVA–Vis reflectance (λ > 0.97 for all body patches) but differed between body regions for NIR reflectance (λ entire wing = 0.98 ± 0.009, λ basal wing = 0.96 ± 0.008, λ thorax = 0.22 ± 0.37).

(b) Evolutionary association between climate and reflectance

NIR reflectance of the thorax and basal wing were positively associated with climate PC1, with higher NIR reflectance in warmer tropical climates and lower NIR reflectance in cooler temperate climates (table 1 and figure 4). For the thorax, however, there was a significant interaction with wing area (table 1 and figure 4b). NIR reflectance increased with increasing PC1 (warmer tropical climates) only for butterflies with smaller wing area. By contrast, UVA–Vis reflectance of the thorax was weakly positively associated with PC1 (no interaction with wing area) and there was no association between UVA–Vis reflectance and PC1 for the basal wing (table 1). The UVA–Vis reflectance of the thorax and basal wing were negatively associated with climate PC2, indicating that the thorax and basal wing area tend to be visibly darker in wetter, cloudier, montane environments (table 1 and figure 4). Notably, the slope and effect size of this relationship were much less than for the relationship between NIR reflectance and climate PC1. Figure 4. Association between environmental variables and reflectivity values for UVA–Vis and NIR regions of the spectra, for two different body regions (thorax: top; basal wing: bottom). Trend line represents the prediction from the model presented in table 1, after accounting for phylogenetic relationships when using the MCC tree. (Online version in colour.)

Table 1. Results of reduced PGLS models predicting UVA–Vis or NIR reflectance across 60 populations (49 species), across 4500 trees. Best models are shown for each body region. NIR models include log UVA–Vis reflectance as a covariate to account for the fact that NIR reflectance is nonlinearly related to UVA–Vis reflectance. The % trees column refers to the percentage of trees where the model presented was the best model. ER refers to the average evidence ratio between the best model (shown) and the null model (for UVA–Vis), and between the best model (shown) and a model with only UVA–Vis reflectance as predictor. Collapse predictors estimate t-value p-value % trees ER UVA–Vis thorax PC1 0.15–0.21 1.27–1.72 0.09–0.20 44 5.03 PC2 −0.36 to −0.28 −2.78 to −2.26 0.005–0.02 basal PC2 −0.37 to −0.30 −2.89 to −2.35 0.004–0.20 99.7 4.09 entire null model 63.4 NIR thorax PC1 1.65–1.78 4.88–5.25 <0.0001 100 >1000 wing size (log) 4.12–4.78 2.70–3.16 0.003–0.008 UVA–Vis (log) 27.98–29.94 7.78–8.39 <0.0001 PC1 × wing −1.74 to −1.58 −3.44 to −3.11 0.001–0.003 basal PC1 0.82–1.03 3.49–4.08 0.0001–0.001 48.3 >1000 wing size (log) 6.95–8.10 2.92–3.57 0.0005–0.004 UVA–Vis (log) 60.5–62.7 15.65–17.08 <0.0001 entire UVA–Vis (log) 46.43–49.55 10.31–11.80 <0.0001 93.2

The significant relationships between either NIR reflectance and PC1 or UVA–Vis reflectance and PC2 were mirrored in the relationship between total reflectance and climate PCs (electronic supplementary material, table S6). Specifically, total reflectance was higher in warmer, tropical climates and lower in cooler, temperate climates, particularly for smaller species (significant interaction between PC1 and wing area; electronic supplementary material, table S6). Total reflectance also decreased in cooler, cloudier climates.

As predicted, in contrast to the thorax and basal wing region, there was no relationship between reflectance of the entire wing and climate PC1 or PC2 (table 1).

There was no association between wing area and climate (PC1: estimate = −0.006–0.002, t-value = −0.41–0.18, p = 0.67–0.99, PC2: estimate = −0.002–0.008, t-value = −0.20–0.69, p = 0.52–0.99).

4. Discussion

The question of whether colour variation across climatic gradients represents an adaptation for thermoregulation remains contentious. Although many studies have shown a relationship between colour and climate, there are multiple possible reasons apart from thermal functions [8]. These include the role of melanin for protection from UV and bacterial damage [48–50], mimicry [51], camouflage [52] and other developmental effects; for example, the polyphenic response to the larval-rearing environment [53,54]. NIR wavelengths are not visually perceived by most terrestrial animals [55], and thermal infrared receptors found in some animals, such as some snakes, beetles, ticks and mites, are largely insensitive to NIR wavelengths [56–59]. Thus, in contrast to reflectance visible to most animals, from about 300 to 700 nm, variation in NIR reflectance should primarily reflect selection for thermoregulation [5]. Our results provide evidence that selection for thermal benefits has shaped the reflectance properties of butterflies. NIR reflectance of the basal wing and thorax, the latter particularly for smaller species, was higher in warmer, tropical climates (lower reflectance in cooler, temperate climates). As predicted, this association was weaker (in terms of both slope and effect size) or absent for UVA–Vis reflectance and reflectance of the entire wing. However, both the thorax and basal wing were also visibly darker in wetter cloudier climates and lighter in drier sunnier climates, supporting previous evidence of climate gradients in melanism in butterflies (e.g. [60–62]), and insects more generally [6,63,64].

Lower reflectance (higher absorption of solar radiation) in cooler, temperate environments may enable butterflies to more rapidly reach high enough temperatures to fly. Flight time strongly limits female reproductive success in high altitude Colias butterflies [62,65] and darker individuals have higher heating rates in a range of butterfly species [16,66–68]. Temperature regulation of the thorax, in particular, is crucial for all flight activities because it is where the flight muscles are located and is the primary site of absorption of incident solar radiation [16]. On the other hand, increased NIR reflectance in warmer and sunnier environments may reduce the chances of over-heating. Heat stress has been shown to reduce female egg production and male lifespan in Colias butterflies, and constrains morphological adaptation, including melanisation, at high altitudes [68–70]. Although butterflies are able to avoid solar radiation by seeking shade or closing their wings to shield their body [67], the optical microstructures on their thorax may nevertheless increase relative NIR reflectance to help reduce the likelihood of exceeding their critical thermal maximum. This is important for all types of basking butterflies, because the thorax and dorsal regions are exposed to the sun during flight activities.

Interestingly, the relationship between NIR reflectance of the thorax and climate was weak or absent for larger species. This was contrary to what we expected because the thermal effects of reflectance tend to increase rather than decrease with body size [9,71]. One explanation may be that the visual relevance of the thorax becomes more important with increasing size (altering NIR properties as a side effect). Another explanation is size-related differences in behaviour [72]. For example, within butterfly communities of tropical woodland–rainforest ecosystems of northeast Australia, darker and larger butterflies tend to be active in the shade and during crepuscular hours, while lighter and smaller butterflies are more active in the sun and midday hours, when NIR reflectivity will have the greatest impact on heat gain [72].

Our results showed that the association between thorax reflectance and climate is similar for the basal wing, which is probably linked to the morphology and physiology of this wing region. There is a low volume of haemolymph fluid in butterfly wings, making them relatively poor conductors of heat [16], and circulation through the entire wing may be too slow to significantly affect thoracic temperature [13–16]. However, the basal region (proximal third of the wing to the thorax) is where most of the veins and haemolymph circulation in the wings are found. Papilio species can increase the temperature of the thorax above ambient air temperature by 40–50% when their wings are fully open and most of this increase can be attributed to the basal wing region [13,14]. However, there is no significant difference in heat gain between living and freshly dead specimens, suggesting that the flow of haemolymph circulation in the basal wing may still not be strong enough for vascular heat transfer, but rather that the radiant energy absorbed by the basal wing is conducted through tissue to the immediately adjacent thorax [12]. This is also influenced by the dense covering of hairs on this region together with dark wing bases characteristic of most dorsally basking butterflies in temperate climates [14].

In contrast to the thorax and basal wing, reflectance of the entire wing was not associated with climate. Wings are the largest and most conspicuous body region, and populations are likely to experience varying selection for optical functions, especially signalling, which may trade off against thermoregulation [73,74]. In our dataset, UVA–Vis (but not NIR) reflectance of the entire wing was higher in males, whereas UVA–Vis and NIR reflectance of the thorax and basal wing were consistently higher in females, suggesting sex differences in selection on reflectance. Furthermore, variation in behaviour may obscure general correlations between wing reflectance and climate, even if wing reflectance enhances thermoregulation. For example, white butterfly species from the Pieridae family alter the angle of their wings to intercept solar radiation and concentrate it onto the thorax [75,76]. Thus, both high reflectance (light coloration) and low reflectance (dark coloration) may enhance heat gain, depending on how wings are used. Experimental data, combined with biophysical models (e.g. [77,78]), are required to assess how entire wing reflectance might influence body temperature under different behavioural and environmental scenarios.

5. Conclusion

By assessing climatic gradients in NIR reflectance, we provide evidence that selection for thermal benefits has shaped the evolution of reflectance properties in Australian butterflies. Our results suggest that butterflies may modulate UVA–Vis and NIR reflectance properties differently to accommodate selection for both optical (e.g. camouflage, signalling) and thermal functions, despite potential constraints imposed by the intercorrelation of wavelengths across the spectrum. Adaptations to control heat gain, such as enhanced NIR reflectivity, may become increasingly important and under stronger selection as the climate warms [63,79]. We also present an example of the use of calibrated NIR digital imaging to study biodiversity and thermal adaptations. Calibrated digital photography is now widely used in ecology and evolutionary biology because it is a relatively inexpensive, accessible and does not require handling of the subject [20,80]. However, its use has been restricted to UVA and visible wavelengths. The application of calibrated NIR imaging may reveal unappreciated variation in reflectance properties and novel thermal adaptations. As butterflies are important models for bioinspired nanophotonic materials [75,81,82], characterizing NIR reflectance properties could inform the development of new functional materials, such as coatings for passive radiative cooling and improved solar energy-capture technology [75,83,84].

Data accessibility

The datasets and trees used for all analyses in the current study are available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.bq862p1 [85].

Authors' contributions

J.T.M., D.S.-F., M.R.K. and K.W. designed and guided the study; J.T.M., A.G.D., J.G. and K.J.R. performed imaging set-up and testing; J.T.M. and K.W. collected all data; M.K. collated the climate data; A.M. reconstructed the phylogeny; J.T.M. and I.M. analysed the data; I.M. and K.J.R. prepared figures; J.T.M. and D.S.-F. wrote the first draft of the manuscript, and all authors contributed substantially to revisions.

Competing interests

We declare we have no competing interests.

Funding

Funding was from a University of Melbourne McCoy Seed Funding grant to D.S.-F. and K.W. D.S.-F. was supported by the Australian Research Council (FT180100216). I.M. was supported by a University of Melbourne McKenzie Fellowship.

Acknowledgements We thank Michelle Dan for providing the illustration in electronic supplementary material, figure S4. We are grateful to John A. Endler and Martin Stevens for access to custom Matlab scripts, Peter Lillywhite for assistance with museum collection access and photography, and Elizabeth Newton for assistance with equipment calibration.

Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.4414733.