Cool habitats support darker and bigger butterflies in Australian tropical forests

Abstract Morphology mediates the relationship between an organism's body temperature and its environment. Dark organisms, for example, tend to absorb heat more quickly than lighter individuals, which could influence their responses to temperature. Therefore, temperature‐related traits such as morphology may affect patterns of species abundance, richness, and community assembly across a broad range of spatial scales. In this study, we examined variation in color lightness and body size within butterfly communities across hot and cool habitats in the tropical woodland–rainforest ecosystems of northeast Queensland, Australia. Using thermal imaging, we documented the absorption of solar radiation relative to color lightness and wingspan and then built a phylogenetic tree based on available sequences to analyze the effects of habitat on these traits within a phylogenetic framework. In general, darker and larger individuals were more prevalent in cool, closed‐canopy rainforests than in immediately adjacent and hotter open woodlands. In addition, darker and larger butterflies preferred to be active in the shade and during crepuscular hours, while lighter and smaller butterflies were more active in the sun and midday hours—a pattern that held after correcting for phylogeny. Our ex situ experiment supported field observations that dark and large butterflies heated up faster than light and small butterflies under standardized environmental conditions. Our results show a thermal consequence of butterfly morphology across habitats and how environmental factors at a microhabitat scale may affect the distribution of species based on these traits. Furthermore, this study highlights how butterfly species might differentially respond to warming based on ecophysiological traits and how thermal refuges might emerge at microclimatic and habitat scales.

In addition to color, large body size also provides physiological benefits for living in cold environments due to low convection rates and higher heat capacities (Porter & Gates, 1969). As a consequence, body size of ectotherms is often negatively correlated with environmental temperature, resulting in body size clines across thermal gradients (Atkinson & Sibly, 1997;Kingsolver & Huey, 2008;Moreno Azócar et al., 2015;Partridge & French, 1996). Because many morphological and physiological traits are linked to climate and the surrounding environment, recent evidence has emerged showing that climate change is triggering increased color lightness and decreased body size for numerous ectotherm species (Angilletta, Niewiarowski, Dunham, Leaché, & Porter, 2004;Gardner, Peters, Kearney, Joseph, & Heinsohn, 2011;Zeuss et al., 2014; but see Connette, Crawford, & Peterman, 2015).
One possible mediator of climate change impacts on biological communities is microhabitats, which are generally decoupled from macroclimatic gradients and thus offer unique thermal regimes for organisms to persist in situ under climate change (Scheffers, Edwards, Diesmos, Williams, & Evans, 2014;. Differences in microclimates within and among habitats, such as closed-canopy forests and open habitats, can be comparable to or greater than gradients across altitude and latitude (Huey et al., 2009;Mark & Ashton, 1992;Scheffers et al., 2013). For small ectotherms that operate at the scale of microhabitats, microscale climate systems are especially important (Bonebrake, Boggs, Stamberger, Deutsch, & Ehrlich, 2014;Pincebourde & Casas, 2015;Potter, Arthur Woods, & Pincebourde, 2013). Another important mechanism in efficient thermoregulation is the adjustment of activity times during the day or throughout the year (Porter, Mitchell, Beckman, & DeWitt, 1973;Stevenson, 1985), which entails finding temperatures at both macro-and microscales that optimize performance (Clusella-Trullas, Blackburn, & Chown, 2011;Grant, 1990;Kearney, Shine, & Porter, 2009).
In this study, we examined the role of temperature and butterfly morphology across tropical forest habitats in north Queensland, Australia, to further understand how thermal gradients at microclimate and habitat scales interact with traits of species assemblages.
They therefore serve as model organisms for testing the consequences of morphology in structuring community traits at micro-and macrohabitat scales and under climate change (Bonebrake et al., 2014;Kingsolver, 1985;Kingsolver & Buckley, 2015). In this study, we first examined species compositions across disparate environments (hot and cool) and linked species distributions in closed and open forest sites to color and body size morphology. We then experimentally examined how color lightness and body size of mounted specimens affect their body temperature change under controlled conditions. Finally, we assessed the influence of evolutionary history on these color traits by constructing a phylogeny and analyzing the effect of phylogenetic and species-specific contributions to measured trait values. The results provide insights into how the environment interacts with morphology to structure communities across microclimatic and habitat scales as well as how that variation could have important implications for how biodiversity will respond to climate change.

| Field surveys
We sampled butterflies in primary rainforests and open woodland habitats of the Australian wet tropical (AWT) bioregion in northeastern Queensland, Australia. Within this region, we conducted our sampling at two locations: Daintree Rainforest National Park (15°57ʹS, 145°24ʹE) and Shiptons Flat (15°42ʹS, 145°13ʹE), from 20 October to 1 Novembr 2014. We used five primary rainforest sites in the Daintree and two primary rainforest sites in Shiptons Flat. We also sampled in two open woodland sites in Shiptons Flat. We used hand nets and binoculars to survey active butterflies in crepuscular (7:00-10:00 a.m. and 3:00-6:00 p.m.) and midday hours (10:00 a.m.-3:00 p.m.) in each habitat (open woodland and primary rainforest). We surveyed along 0.5-km transects for 30-minute intervals, and we marked each captured butterfly to avoid double counting during each sampling. For those individuals not easily identified during sampling, we collected the specimen for identification in the laboratory. Additionally, we collected and mounted at least one specimen for each species for color lightness analysis.
We recorded the time of each capture and categorized the spot of capture as either sunny or shady. We used one Thermochron iButton data logger (model: DS1921) in each habitat to record the ambient temperature throughout the sampling period. Each temperature logger was suspended approximately 1.5 m above the ground and hung beneath a plastic funnel to avoid direct sunlight (Scheffers et al., 2013;Shoo, Storlie, Williams, & Williams, 2010).

| Collection of morphology information
We photographed specimens using standardized settings (exposure time: 1/200 s, ISO speed: 100, aperture: F/16) and a flash from a fixed distance to record color. Basking postures for absorbing heat may vary by species, so we took photographs of both the dorsal and ventral sides of the wings for each specimen. The lightness value of each species was obtained using Adobe Photoshop CC 2014 software.
Because studies have shown that the basal part of the wing and body are important for thermoregulation (Kingsolver, 1985;Wasserthal, 1975), we analyzed the color of the body, together with the one-third wing area immediately adjacent to the body. We also analyzed the color of the body and the whole specimen (whole wing area plus body) to test the consistency of the lightness value across regions of the specimens. We used the quick selection tool in Photoshop to choose each region for color analysis and averaged the chosen region using the average blur filter. We used the mean of the value of the red, green, and blue channels as the final estimated lightness value (a value between 0 and 255-low values denote dark butterflies and high values indicate lighter individuals). For most of the species, we picked one best-preserved specimen to represent the morphology. For species with multiple specimens, we took the average color lightness of all specimens, averaging male and female lightness values when specimens of both sexes were available. As wingspan is known to be an effective and convenient proxy for body size (Sekar, 2012), we collated wingspan (the distance between the tips of the forewings) data from the Butterflies of Australia field guide as an indicator of body size (Braby, 2000).

| Thermal experiment
Using the best-preserved collected mounted specimens for 26 species, we conducted an experiment with a FLIR thermal infrared camera (model: E6) over the course of two sunny days from 12 to 3 p.m.
Our E6 model records 19,200 samples of temperature per photograph (160 × 120 pixels). This time period was chosen to ensure stable and strong solar radiation. We standardized our start temperatures by placing mounted specimens in full shade and exposed each specimen to full sun for 2 min. Using the thermal camera, we recorded the surface temperature of the thorax every 30 s. In this way, we collected starting ambient shade temperature of the thorax in addition to four sun-exposed thorax temperatures.
We kept all other priors as default except cold Markov chain with a temperature parameter of 0.16, and we adjusted the mean branch length prior to 0.01 (brlenspr = unconstrained:exponential(100.0)) to reduce the likelihood of stochastic entrapment in local tree length optima (Brown, Hedtke, Lemmon, & Lemmon, 2010;Marshall, 2010).
The resulting standard deviation of split frequencies was <0.01. The parameters sampled during the MCMC were imported into Tracer v.1.5 showing adequate effective sample size (>200, Rambaut & Drummond, 2007). We excluded 4,000 initial samples (around 33%) as burn-in of each MCMC run from the summary analysis, and we calculated a 50% majority-rule consensus tree from the post-burn-in trees. The posterior probabilities (PPs) of each node were summarized (Larget & Simon, 1999), and we used these to infer support for individual clades: nodes with PP values ≥.95 were regarded as well or strongly supported (Yang & Rannala, 1997).

| Phylogenetic analysis
We extracted 1,000 trees randomly from the post-burn-in trees in one of the MCMC runs. We removed the 16 taxa which were added for phylogenetic reconstruction ( (Maddison & Maddison, 2011). For each sampled tree, we used Pagel's λ (Pagel, 1999) to calculate the phylogenetic signal of traits with the phylosig function in the R package phytools (Revell, 2012). We then rescaled the traits value by Pagel's λ using the rescale function in the R package geiger (Harmon, Weir, Brock, Glor, & Challenger, 2008) and used Lynch's comparative method (Lynch, 1991) to separate the total trait values to phylogenetic components and specific components for the 1,000 trees with function compar.lynch in R package ape (Paradis, Claude, & Strimmer, 2004). The phylogenetic components represent the ancestral contribution to the trait, while the specific components represent the species-specific variance of the trait (Zeuss et al., 2014).
In this way, 1,000 results for phylogenetic components and specific components for each species were generated.

| Community analysis
For all species for which we obtained abundance and environmental information, we used individual-based rarefaction curves to compare the species richness of rainforest and woodland habitats. We also conducted a nonmetric multidimensional scaling (NMDS) to examine community-level differences across habitats and sites, using each sampled location as a replicate. We excluded one rainforest site in the Daintree (Mt. Sorrow) from the community analysis due to the extreme low abundance of butterflies. We then conducted a PERMANOVA test using the adnois function in the R package vegan (Oksanen et al., 2013) to quantify the difference in community composition across habitats.

| Temperature gradient
We extracted the ambient temperature for each point in time at which a butterfly was observed or captured using data from the iButton data loggers. We then applied a multiple linear regression model to estimate the temperature gradient experienced by butterflies across habitat and time (crepuscular and midday hours) with the predict function in the R package stats (R Core Team, 2013).

| Thermal experiment
We applied linear mixed-effects models to estimate the contribution of color lightness and body size to the observed increase in body temperature during the thermal experiment. As the thermal experiment was conducted over 2 days, we used day as a random factor and color lightness and wingspan as fixed factors. We used the lme function from the R package nlme (R Development Core Team, 2013).
We performed a linear correlation of color lightness between whole wing with body, basal wing with body, and body only to test the consistency of color lightness across these body regions. We constructed two models with and without an interaction between color lightness and wing size and used Akaike information criterion (AIC) model selection to choose the model with the lowest AIC values for analyses (Akaike, 1974).

| Relationship between butterfly morphology and environmental factors
To relate habitat and assemblage trait characteristics, for each individual, we used the species-specific morphological data as well as the environmental data (habitat, sun vs. shade) recorded at the time when the individual was observed. We applied a multiple linear regression to analyze how butterfly color lightness and body size was affected by habitat, time, and sunny/shady conditions. We also constructed multiple linear regression models with and without interactions between different environmental factors and applied AIC to choose the model with the lowest AIC value as the best-fit model for morphology and environmental gradient analyses. We then used the predict function in the R package stats (R Development Core Team, 2013) to estimate butterfly color lightness and body size across the gradients generated by different factors based on the best selected models. We used the whole dataset (46 species with 408 individuals) for the analysis above.
For the 1,000 phylogenetic components and specific components generated by the 1,000 phylogenetic trees for traits of each species (see the section 2.5, phylogenetic analysis), we separately ran the same selected multiple linear regression model for each component against environmental factors. We then used the predict function to estimate phylogenetic and specific components at each of the gradients based on the linear regression model results and averaged across the 1,000 estimates. Due to a lack of phylogenetic information for some rare species, we used a subset of the dataset (37 species with 363 individuals) for the phylogenetic analysis. All analyses were conducted in R v.3.2.2 (R Development Core Team, 2013).
Although our focus in this analysis was on individuals, we also used species as a unit to look into the presence-absence of different species traits across environmental gradients by applying the lowest AIC model for traits and environmental factors. Finally, in order to assess how variation in the abundance of common species across gradients might be affecting the observed patterns, we chose the top five mostabundant species separately from our rainforest and woodland habitats and plotted the abundance of those species at each combination of light exposure, time of day, and habitat together with corresponding morphological data.

| Summary of butterfly assemblages
We recorded 408 individuals of 46 species through capture and binocular observation (Figs S1 and S2). Species richness was higher in rainforest (38 species) than in woodland (19 species; 11 species were found in both habitats; Fig. S3). The NMDS plot demonstrated some differentiation of the community composition between rainforest and woodland but no significant difference by the PERMANOVA test (F = 1.15, p = .17; Fig. S4).
The phylogenetic relationships of the species were summarized in Fig. S5. The phylogenetic relationship between different families in our study was similar to previous work (Heikkilä, Kaila, Mutanen, Peña, & Wahlberg, 2012;Regier et al., 2013). Although the relationship of the out-group was ambiguous, the "true butterflies" were shown to be monophyletic (PP = .98). All deep nodes were strongly supported. The families Nymphalidae and Hesperiidae form a strong supporting clade sister to the clade containing the families Pieridae and Lycaenidae.
Papilionidae is basal to the rest of the families.

| Temperature gradients
Ambient temperature experienced by butterflies was significantly different between habitats and time of sampling (R 2 = .68, Fig. S6).
The temperatures experienced in closed rainforests were approximately 2.3 ± 0.3°C (mean ± SE) lower than those in open woodlands, and temperatures experienced during crepuscular hours were approximately 3.4 ± 0.2°C (mean ± SE) lower than those experienced during midday hours. The coolest temperature observed in our study was approximately 25°C during crepuscular survey hours in rainforests, while the hottest temperature recorded was approximately 33°C during midday survey hours in open woodlands (Fig. S6).

| Thermal experiment
Color lightness of three regions of the body and wing was highly correlated as indicated by linear regressions (p < .001, R 2 > .8). Given that the basal wing and body are considered the most important parts for butterfly thermoregulation (Kingsolver, 1985;Wasserthal, 1975), we focused on this region for subsequent analyses and note that here and elsewhere in the results and discussion "color lightness" refers to that of the basal part of the wing plus body. Our most parsimonious linear model suggests that color lightness negatively affected the rate of butterfly body warming under exposure to solar radiation, while wingspan positively affected warming (Table 1).

| Butterfly color lightness and body size across habitat gradients
Based on the model with the lowest AIC value (Table 2), we found that for both dorsal and ventral sides of the wings, color lightness of butterfly individuals predicted in closed rainforest was significantly darker than that found in open woodlands (Fig. 1) (Table 3). Body sizes of butterfly individuals in closed rainforests were also predicted to be larger relative to those in open woodlands ( Fig. 2 and Table 3).
Within habitat, we observed a temporal and spatial effect on butterfly lightness. Specifically, in rainforests, individuals active during midday hours were lighter than those active during crepuscular hours ( Fig. 1 and Table 3). In addition, individuals active in the sun were lighter than those in the shade, and this pattern became more pronounced during hotter midday conditions (Fig. 1). However, these patterns were not apparent in the woodland (Fig. 1). We found no significant effect of activity time and sun/shade level on wingspan ( Fig. 2 and Table 3).
We found no significant difference in color lightness and size between microhabitats when data were analyzed at the species level (Figs S7 and S8). However, the abundance of the most-abundant species for each habitat (five species from each habitat with only one abundant in both habitats: Mycalesis terminus), which account for near 70% of total abundance, varied by habitat and time of day (Figs 3 and   4). The most-abundant species with darker colors are abundant under cool conditions such as in closed rainforests and during crepuscular hours (Fig. 3). Meanwhile, the most-abundant species with relatively larger body sizes are abundant in closed rainforests and conversely those with relatively smaller body sizes are abundant in hot conditions such as open woodland and midday hours (Fig. 4).
Wingspan and color lightness were both highly phylogenetically correlated as lambda was greater than 0.9 for both traits (Fig. S9).
Analyses of both phylogenetic components and specific components of color lightness and wingspan of butterfly individuals revealed that darker coloration and larger sizes were predicted in closed rainforest (Figs 5 and 6). Moreover, analyses of color lightness components showed that individuals active in the sun and midday were lighter than those active in the shade and crepuscular within rainforest. The pattern for specific components of color lightness across environment gradient was stronger for butterfly ventral side than for the dorsal side (Fig. 5).

| DISCUSSION
Our results suggest that darker and larger butterfly individuals tend to prefer cooler conditions both within and across habitats in tropical Australia. We observed consistent patterns of lightness and size across each level of exposure in space and time, spanning microhabitat (shade to sun), time of day (crepuscular to midday), and habitat (rainforest to open woodland). Although these patterns are highly influenced by phylogeny, the strong correlation between environmental factors and specific components shows that temperature gradients across habitats play an important role in structuring butterfly T A B L E 1 Multiple linear regression models of color lightness and body size effects (mean ± SE) on experimental body temperature increase (lightness = color lightness of basal wing and body, size = wingspan)   traits. The increased prevalence of dark and large individuals in cool conditions is supported by the results from our thermal experiment, which showed that dark bodies and larger sizes can accelerate heat gain, a finding consistent with other published research (Kingsolver, 1985;Porter & Gates, 1969;Wasserthal, 1975;Zeuss et al., 2014).
Our findings support two prominent ecological hypotheses: that melanism scales with temperature (thermal melanism hypothesis) and that body size scales with temperature (Bergmann's rule) (Bergmann, 1847;Clusella-Trullas et al., 2007;Gardner et al., 2011;Moreno Azócar et al., 2015;Partridge & French, 1996). Interestingly, we validated these two patterns at the interspecific level in the field-a finding which complements prior field research which has largely been focused on intraspecific variances at the population level, and macroecological research which has focused on large-scaled distribution databases (Davis, Chi, Bradley, & Altizer, 2012;Ellers & Boggs, 2004;Moreno Azócar et al., 2015).
Temperature clearly differed between rainforest and open woodland sites even though these two habitats share a common ecotone. Our results suggest that trait variation can be shaped by these small-scale thermal regimes (Baudier, Mudd, Erickson, & O'Donnell, 2015;Duffy, Coetzee, Janion-Scheepers, & Chown, 2015)-a finding that is consistent with recent studies at the microhabitat (Baudier et al., 2015;Kaspari, Clay, Lucas, Yanoviak, & Kay, 2015) and ecosystem scale (e.g., temperate ecosystems, Moreno Azócar et al., 2015;Zeuss et al., 2014). Studies that focus solely on macroscale environmental gradients may be insufficient in resolution to detect ecotypic trait patterns of small-bodied ectotherms, especially in complex environments such as tropical forests (Duffy et al., 2015). Solar radiation also likely plays a significant role in shaping the pattern of butterfly color lightness across and within habitats. For example, within rainforests, we observed clear differences in predicted color lightness between butterfly individuals in sun versus shade microhabitats as well as large differences in lightness between crepuscular and midday hours when solar radiation was strongest. In open woodlands, as expected, we did not observe any obvious difference in predicted color lightness of butterfly individuals between sun/shade conditions or crepuscular/midday hours, which is likely due to the overall uniformity of heat and radiation exposure in these open habitats. Therefore,  relative to the closed rainforest, the comparatively homogenous environment of open woodlands yields lower variability in traits such as color lightness. Across the tropics, large areas of rainforests continue to be selectively logged (Hansen, Stehman, & Potapov, 2010;Peh et al., 2011;Pert, Butler, Bruce, & Metcalfe, 2012). This not only undermines critical life-cycle processes, but also exposes these communities with long evolutionary histories of low-light environments to high levels of solar radiation and temperature, for which they are neither acclimated nor adapted to. Our findings suggest that deforestation and selective logging could reduce the availability of thermally buffered microhabitats for tropical species (Gardner et al., 2009;Huey et al., 2009; and dramatically increase their exposure to direct solar radiation (Brown, 1993;Carlson & Groot, 1997). Both human and natural disturbances can open the canopies of rainforests, which can present physiological challenges to closed-forest butterflies (Koh, 2007).
Morphological traits are the result of "long-term" evolutionary influences (phylogenetic component) as well as more recent adaptations or responses to the environment (specific component) (Zeuss et al., 2014). The effect of historic climate on biological patterns in the AWT is well understood (Schneider, Williams, Bermingham, Dick, & Moritz, 2005;Williams, Bolitho, & Fox, 2003) with all analyses conducted to date showing that fluctuations in rainforest extent during the Quaternary Period has been the single largest determinant of current regional-scale patterns of vertebrate assemblage structure and biodiversity (Graham, VanDerWal, Phillips, Moritz, & Williams, 2010;Schneider, Cunningham, & Moritz, 1998;Schneider & Moritz, 1999;Schneider et al., 2005). Importantly, the predominant effect of midlate Pleistocene climate fluctuations was extinction rather than speciation with many communities effectively filtered by advantageous or disadvantageous traits (Moritz, Patton, Schneider, & Smith, 2000;Williams & Pearson, 1997).
We partitioned the variance in color lightness and wingspan into phylogenetic and species-specific components in order to determine the influence of phylogenetic inertia on our observed patterns of color lightness and body size and found that both the phylogenetic and species-specific components showed differences in traits between hot and cold time periods and spatial scales. The phylogenetic component is manifested as different species prevalence between rainforest and woodland habitats, whereas the specific components of the mea- convincing support that human-induced disturbances that change ambient temperature and exposure to solar radiation could filter communities by these traits (Zeuss et al., 2014). Tropical ectotherms have high projected risk of exceeding thermal limits under climate warming (Bonebrake & Deutsch, 2012;Huey et al., 2009;Tewksbury, Huey, & Deutsch, 2008), and traits that are beneficial under cool environments may lose their value or become costly as the climate warms (Willis et al., 2015). Behavioral studies suggest that butterflies need to keep their body temperatures within a range of 30-40°C for flight (Kingsolver, 1983;Kingsolver & Watt, 1983;Watt, 1968 (Bonebrake et al., 2014).
There are other factors that structure butterfly traits at both local and broad scales which were not considered in our study but should be acknowledged to facilitate complementary research. For example, in this study, we focused on morphology at the species and community/ assemblage level rather than on variation among individuals and sexes (e.g., Ellers & Boggs, 2004). A more detailed trait analysis considering intraspecific variance in lightness and size could also provide valuable information on how plastic these traits might be in variable environments. Moreover, the color lightness in this study only captures variation in pigment coloration, yet structural or iridescence might also influence thermoregulation and fitness in butterflies (Tamáska et al., 2013). Additionally, biotic factors such as predation and competition exert selective pressures on species morphology and can shape the morphological pattern of assemblages across environmental gradients (Chai & Srygley, 1990;Cook, Mani, & Varley, 1986). This is es- Our study emphasizes the important role of microclimates and habitat in shaping the color lightness and body sizes of tropical butterflies. Within a phylogenetic framework, we found bigger and darker butterflies in cool microhabitats, suggesting that these traits confer some benefit for ectotherms in cool environments at fine spatial scales, even in lowland tropical environments. These results together highlight the complex interactions between morphology and microclimate, across ecological and evolutionary time scales, and how species vulnerability to climate change could be mediated by these interactions.

ACKNOWLEDGMENTS
We thank our field research assistant Ms. Nadiah Roslan and all the volunteers who have helped in data collection. Prof. Richard Saunders and his laboratory contributed generously in the analysis and interpretation of the phylogenetic data. We also thank Mr. Toby Tsang and two anonymous reviewers for their valuable comments on the manuscript.

FUNDING INFORMATION
This work was generously supported by the Centre for Tropical