Altitude and life‐history shape the evolution of Heliconius wings

Abstract Phenotypic divergence between closely related species has long interested biologists. Taxa that inhabit a range of environments and have diverse natural histories can help understand how selection drives phenotypic divergence. In butterflies, wing color patterns have been extensively studied but diversity in wing shape and size is less well understood. Here, we assess the relative importance of phylogenetic relatedness, natural history, and habitat on shaping wing morphology in a large dataset of over 3500 individuals, representing 13 Heliconius species from across the Neotropics. We find that both larval and adult behavioral ecology correlate with patterns of wing sexual dimorphism and adult size. Species with solitary larvae have larger adult males, in contrast to gregarious Heliconius species, and indeed most Lepidoptera, where females are larger. Species in the pupal‐mating clade are smaller than those in the adult‐mating clade. Interestingly, we find that high‐altitude species tend to have rounder wings and, in one of the two major Heliconius clades, are also bigger than their lowland relatives. Furthermore, within two widespread species, we find that high‐altitude populations also have rounder wings. Thus, we reveal novel adaptive wing morphological divergence among Heliconius species beyond that imposed by natural selection on aposematic wing coloration.

Identifying the selective forces driving phenotypic divergence among closely related species lies at the core of evolutionary biology research. Adaptive radiations, in which descendants from a common ancestor rapidly fill a variety of niches, are ideal systems to investigate morphological divergence (Schluter 2000). The study of adaptive radiations has revealed that evolution often comes up with similar solutions for similar problems at the phenotypic and genetic levels (Losos 2010;Marques et al. 2019). Speciose groups that have repeatedly and independently evolved convergent adaptations to life-history strategies and environments are good systems in which study selection drivers (Schluter 2000). Nevertheless, adaptive phenotypic evolution is often complex and multifaceted, with more than a single selective force in action (Maia et al. 2016;Nosil et al. 2018). For example in birds, sex differences in plumage coloration are driven by intraspecific sexual selection, while natural selection drives sexes toward more similar colorations (Dunn et al. 2015). Integrative approaches that make use of tractable traits across well-resolved phylogenies are needed to explore the selective forces driving phenotypic evolution.
Butterfly wing coloration has been the focus of considerable research effort and major strides have been made toward understanding how and when evolution leads to complex wing color patterns, conferring aposematism, camouflage, or a mating advantage (Merrill et al. 2012;Chazot et al. 2016;Nadeau et al. 2016). The dazzling diversity of butterfly color patterns among species has perhaps obscured the less conspicuous phenotypic diversity of wing shapes and sizes, which are more often regarded as the result of sexual selection, flight trade-offs or developmental constraints (Singer 1982;Allen et al. 2011), rather than drivers of local adaptation and species diversification (Srygley 2004a;Cespedes et al. 2015;Chazot et al. 2016). A recent review assessing the ecology of butterfly flight, identified habitat, predators, and sex-specific behaviors as the selection forces most likely driving wing morphology variation, but highlighted the need for further phylogenetic comparative studies that identify the adaptive mechanisms shaping wings (Le Roy et al. 2019).
Differences in behavior between sexes have been identified as one of the main drivers of wing aspect ratio and size sexual dimorphism in insects (Rossato et al. 2018a;Le Roy et al. 2019). In butterflies, males tend to spend more time looking for mates and patrolling territories, while females focus their energy on searching for suitable host plants for oviposition (Rossato et al. 2018b). The same wing trait can be associated with different life history traits in each sex, resulting in sex-specific selection pressures. For example, in the Nearctic butterfly Melitaea cinxia, wing aspect ratio only correlates with dispersal in females, as males experience additional selection pressures that counteract selection for dispersal wing phenotypes (Breuker et al. 2007). Sex-specific behaviors can impact wing aspect ratio and size, but differences in life histories, even across closely related species, could also have large impacts on the strength and direction of these effects (Cespedes et al. 2015;Chazot et al. 2016).
Another important source of phenotypic variation in insect wings is the physical environment they inhabit throughout their range. Air pressure decreases with altitude, which in turn reduces lift forces required for flight. To compensate for this, insects may increase wing area relative to body size to reduce the velocity necessary to sustain flight (Dudley 2002;Dillon et al. 2018). Wing aspect ratio in Drosophila melanogaster has been observed to vary adaptively across latitudes and altitudes, with wings getting rounder and larger in montane habitats, possibly to maintain flight function in lower air pressures (Stalker and Carson 1948;Pitchers et al. 2012;Klepsatel et al. 2014).
In butterflies, high aspect ratios, that is, long and narrow wings, reduce drag caused by wing tip vortices, thus lowering the energy required for flight and promoting gliding for longer distances (Le Roy et al. 2019). Variation in wing phenotypes can occur at the microhabitat level, for example, Morpho butterfly clades in the understory have rounder wings than canopyspecialist clades, presumably for increased maneuverability (Chazot et al. 2016). An extreme case of environmental effects on wing morphology can be found in Lepidoptera inhabiting the windy, barren highlands of the Andes, where an interaction between behavioral sex differences and extreme climatic conditions have led to flightlessness in females of several species (Pyrcz et al. 2004).
Heliconius is a genus of Neotropical butterflies that has been studied for over two centuries with a well resolved phylogeny (Kozak et al. , 2018. It represents a striking case of Müllerian mimicry, with co-occurring subspecies sharing warn-ing wing color patterns to avoid predators and leading to multispecies mimicry rings across South America (Merrill et al. 2015). Wing aspect ratio and size are part of the mimetic signal (Jones et al. 2013;Mérot et al. 2016;Rossato et al. 2018a). Wing morphology is involved in many aspects of Heliconius biology other than mimicry, such as mating or flight mode, but these have been less well studied (Rodrigues and Moreira 2004;Srygley 2004a;Mendoza-Cuenca and MacÍas-Ordóñez 2010). As the only butterflies that pollen-feed, their long life-spans and enlarged brains allow them to memorize foraging transects that are repeated daily following a short dispersal post-emergence phase of up to 1.5 km (Cook et al. 1976;Jiggins 2016).
Larval gregariousness has evolved independently three times across the phylogeny, with some species laying clutches of up to 200 eggs, while others lay eggs singly and larvae are often cannibalistic (Beltrán et al. 2007). Gregarious Heliconius species would be predicted to have larger sized females to carry the enlarged egg load, as is the case with most Lepidoptera (Allen et al. 2011). Another striking life history trait is pupal-mating, which is only found in one of the two major clades (hereafter the "erato clade"), having arisen following the most basal split in the Heliconius phylogeny. This mating strategy involves males copulating with females as they emerge from the pupal case (Deinert et al. 1994;Beltrán et al. 2007). Pupal-mating leads to a whole suite of distinct selection pressures but these are hard to tease apart from the effects of phylogeny due to its single origin (Beltrán et al. 2007;Thurman et al. 2018). Further ecological differences could arise from adaptation to altitude. Some species are relatively high-altitude specialists, such as H. telesiphe and H. hierax found above 1000 m, while others range widely, such as H. melpomene and H. erato, which can be found from 0 to 1800 m above sea level (Rosser et al. 2015;Jiggins 2016). Potential adaptations to altitude are yet to be explored.
The wide range of environments that Heliconius species inhabit, together with their diverse natural history and well-resolved phylogeny, make them a good study system for teasing apart the selective forces driving wing phenotype (Merrill et al. 2015;Jiggins 2016). Here, we examine variation in wing aspect ratio and size across 13 species that span most of the geographical range of the Heliconius genus. First, we photographed thousands of wings collected by many Heliconius researchers since the 1990s from wild populations across South and Central America, covering a 2100 m elevation range (Fig. 1A). Wing dimensions for 3515 individuals, obtained with an automated pipeline and standardized images, were then used to address the following questions. (1) Are there size and aspect ratio sexual dimorphisms, and if so, do they correlate with known life-history traits? (2) To what extent are wing aspect ratio and size variation explained by shared ancestry? (3) Are wing aspect ratio and size affected by the elevations species inhabit?

STUDY COLLECTION
The wild specimens studied here were collected using hand nets between 1998 and 2018 in 313 localities across Panama, Colombia, Ecuador, French Guiana, Suriname, and Peru (Fig. 1A), and stored in the Department of Zoology, University of Cambridge (Earthcape database). Collection altitudes ranged from sea level to 2100 m above sea level (Fig. 1A). Detached wings were photographed dorsally and ventrally with a DSLR camera with a 100 mm macro lens in standardized conditions. 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).

WING MEASUREMENTS
Damage to wings was manually scored in all the images and damaged specimens were excluded from our analyses. To obtain wing measurements from the images, we developed custom scripts for Fiji (Schindelin et al., 2012), to automatically crop, extract the right or left forewing, and perform particle size analysis (Fig. 1B). Butterflies predominantly use their forewings for flight (Wootton 1992;Le Roy et al. 2019) and hindwings tend to be more damaged in Heliconius due to in-flight predation and fragile structure, thus we only include forewings here. Forewing and hindwing areas are tightly correlated in this genus (Strauss, 1990). For wing area, we obtained total wing area (in mm 2 , hereafter "size").
For examining wing aspect ratio, the custom scripts first fitted an ellipse to the forewings and measured the length of the longest axis and the length of the axis at 90°to the former (Fig. 1C). Aspect ratio corresponds to the length of the major axis divided by the length of the minor axis, hereafter "aspect ratio" (Fig. 1C). Aspect ratio was used as a proxy of wing shape, as it has been widely used in previous studies and correlates to gliding efficiency (Le Roy et al. 2019). The data were checked for visual outliers on scatter-plots, which were examined, and removed from the analyses if the wing extraction pipeline had failed.

STATISTICAL ANALYSES
All analyses were run in R version 2.13 (R Development Core Team 2011) and graphics were generated with the package ggplot2 (Ginestet 2011). Packages are specified below. All R scripts can be found in the public repository Zenodo (Zenodo: https://doi.org/10.5281/zenodo.3491029), including custom Fiji scripts for wing image analysis. Species and sexes mean trait values were calculated for the 13 Heliconius species in our study. Each species had more than 30 individuals and all individuals had accurate locality and altitude data (Table S1), resulting in a dataset of 3515 individuals.

Sexual dimorphism across species
Sexual dimorphism in wing area and aspect ratio was estimated as the female increase in mean wing area and aspect ratio with respect to males, thus negative values represent larger trait values in males, while positive values represent larger trait values in females. Pairwise t-tests were used to estimate the significance of sexual size/shape dimorphism in each species.
We modeled variation in wing area and aspect ratio sexual dimorphism across species with ordinary least squares (OLS) linear regressions, implemented in the 'lm' function. For models of sexual wing area and aspect ratio sexual dimorphism, predictor variables initially included larval gregariousness of the species (gregarious or solitary, as classified in Beltrán et al. 2007), mating strategy (pupal-mating vs. adult-mating clade), species mean wing aspect ratio and area, and species wing aspect ratio or size sexual dimorphism (respectively). Wing size sexual dimorphism had a marginally significant phylogenetic signal (Abouheif Cmean = 0.25, P = 0.05), so we present the sexual size dimorphism model incorporating phylogeny as correlation term in the Supporting Information (Tables S3 and S4). We used backward selection with Akaike Information Criterion corrected for small sample sizes (AICc, Hurvich and Tsai 1989) where the best models had the lowest AICc values, implemented with the package MuMin (Bartón 2019). We report the overall variation explained by the fitted linear models (R 2 ) and the relative contributions of each explanatory variable (partial R 2 ), estimated with the package relaimpo (Grömping 2006).

Variation across species
To test whether variation in wing aspect ratio and area across species was constrained by shared ancestry, we calculated the phylogenetic signal index Abouheif's Cmean (Abouheif 1999), which is an autocorrelation metric suitable for datasets with a relatively low number of species and that does not infer an underlying evolutionary model (Münkemüller et al. 2012). Observed and expected distribution plots for phylogenetic signal estimates are shown in the Supporting Information and were computed with the package adephylo (Jombart and Dray 2010). We used a pruned tree with the 13 species under study from the most recent molecular Heliconius phylogeny . We plotted centered trait means across the phylogeny with the function barplot.phylo4d() from the package phylosignal (Keck et al. 2016). To test and visualize phylogenetic signal further, we built phylocorrelograms for each trait with the function phyloCorrelogram() of the same package, which estimates Moran's I autocorrelation across matrices with varying phylogenetic weights. Then, the degree of correlation (Morans' I) in species trait values can be assessed as phylogenetic distance increases (Keck et al. 2016).
To study variation in wing area and aspect ratio across species, we took a phylogenetic comparative approach. These methods assume that species-specific mean trait values are a good representation of the true trait values of the species under study; in other words, the within-species variation is negligible compared to the across-species variation (Garamszegi 2014). To test this, we first used an ANOVA approach, with species as a factor explaining the variation of mean trait values. We then estimated within-species trait repeatability, or intraclass correlation coefficient (ICC), with a linear mixed model approach. This requires the grouping factor to be specified as a random effect, in this case species, with a Gaussian distribution and 1000 parametric bootstraps to quantify uncertainty, implemented with the function rptGaussian() in rptR package (Stoffel et al. 2017). By specifying species as a random effect, the latter approach estimates the proportion of total trait variance accounted for by differences between species. A trait with high repeatability indicates that species-specific trait means are reliable estimates for further analyses (Stoffel et al. 2017). We, nevertheless, accounted for within-species variation in the models described below.
To test the effect of altitude on wing aspect ratio and size across species, we used a phylogenetic generalized least squares (PGLS) approach. Species wing trait means may be correlated due to shared ancestry (Freckleton et al. 2002;Chazot et al. 2016). Therefore, to explore the effects of the environment on the traits under study, models that incorporate expected correlation between species are required, such as PGLS. Although often ignored, these models assume the presence of phylogenetic signal on the model residuals of the trait under study (here wing aspect ratio or size) controlling for covariates that affect the trait mean (allometry, sex ratio), and not just phylogenetic signal on the species mean trait values (Revell 2010;Garamszegi 2014). Thus, to check if this assumption was met, we estimated phylogenetic signal as described above (Keck et al. 2016) for the residuals of a generalized least squares (GLS) of models that had wing aspect ratio or size as response variables, and the size and aspect ratio (respectively) and sex ratio as explanatory variables, to ensure this assumption of PGLS model was met. To visually inspect phylogenetic sinal on the residuals, we obtained phylogenetic correlograms for these and centered trait residuals for plotting across the phylogeny as detailed above for trait means (presented in the Figs. S3 and S4; Keck et al. 2016).
Significant phylogenetic signal was detected in mean wing size and in the residuals of both traits, wing aspect ratio, and area regression models ( Figs. S4 and S5), so we used maximum loglikelihood PGLS regression models with the phylogenetic correlation fitted as a correlation term, implemented with the gls() function from the nmle package (Pinheiro et al. 2007). We assumed a Brownian motion model of trait evolution for both traits, by which variation across species accumulated along all the branches at a rate proportional to the length of the branches (Freckleton et al. 2002). To select the most supported model given the available data, that is, one that improves model fit while penalizing complexity, we used the Aikaike Information Criteria corrected for small sample sizes (AICc, Hurvich and Tsai 1989), where the best models had the lowest AICc values, implemented with the package MuMin (Bartón 2019). Maximal PGLS models included species mean altitude and distance from the Equator (to control for potential latitudinal clines), sex ratio in our samples interacting with either wing aspect ratio or wing size, to control for potential allometric and sexual dimorphism relationships, which could be different among closely related taxa (Outomuro and Johansson 2017). Most species are found in the Andean mountains or the Amazonian region near the Equator, so we did not have much power to examine variation with latitude in wing aspect ratio and size across species, but we included distance from the Equator as an explanatory variable in the PGLS models to account for it. Minimal PGLS models consisted of the trait under study explained solely by its intercept, without any fixed effects. All model selection tables can be found in the Supporting Information (Tables S3 and S5). Finally, we weighted PGLS regressions to account for unequal trait variances and unbalanced sample sizes across species (for sample sizes and standard errors of species' trait means see Table S1). This was achieved by modifying the error structure of the model with combined variances obtained with the function varFixed() and specified with the argument "weights" (Pinheiro et al. 2007;Paradis 2012;Garamszegi 2014). In this study, 74.8% of the individuals were collected in the last 10 years, thus we did not have power to detect any changes in wing morphology across species potentially incurred by climate change (Fig. S1). Future studies could focus on temporal changes in wing morphology in areas and species that have been well sampled throughout the years.

Variation within species
We selected the two most abundant and geographically widespread species within our dataset, H. erato (n = 1685) and H. melpomene (n = 912; Table S1), to examine variation in wing area and aspect ratio within species. We modeled variation in size and aspect ratio with ordinary least squares (OLS) linear regressions for each species, implemented in the 'lm' function. For all models, predictor variables initially included the terms altitude, distance from the Equator, longitude, aspect ratio or wing area, and sex, as well as the plausible interactions between them (Table S5). We then used step backward and forward selection based on AIC with the function stepAIC(), from the MASS package (Ripley, 2011;Zhang, 2016; full models and model selection tables in Tables S5 and S6).

Results
We obtained intact-wing measurements for 3515 individuals of 13 Heliconius species from across the phylogeny and from over 350 localities ( Fig. 1; Table S1). We have made all of these wing images publicly available at the Zenodo repository.

SEXUAL DIMORPHISM
Sexual dimorphism in wing area was found throughout the phylogeny, but in opposing directions in different species (Fig. 2).
Mean sizes were significantly or marginally significantly different among sexes in nine species, all of which were represented by more than 40 individuals (Table S2 for two sample T-test summary statistics), indicating that the nonsignificant trends in other species probably reflect a lack of power caused by low numbers of females typically collected in the wild (Table S1). The six species with trends toward larger females have gregarious larvae (pink, Fig. 2), whereas the seven species with trends toward larger males lay eggs individually (black, Fig. 2). Larval gregariousness alone explained 69% of the total natural variation in sexual size dimorphism across species (Table 1; Gaussian LM: F 1,11 = 27.2, P < 0.001, R 2 = 0.69). There was a marginally significant phylogenetic signal in sexual size dimorphism (Abouheif's Cmean = 0.24, P = 0.05; Fig. S3), so we repeated the analysis accounting for phylogeny and the results are presented in the Supporting Information. This would be expected from the evolutionary history of gregariousness, as it is present in all species of three lineages that are well represented in our study (Beltrán et al. 2007). However, when accounting for phylogenetic correlation in the model larval gregariousness remained a significant predictor of size sexual dimorphism (Table S4).

PHYLOGENETIC SIGNAL
The 13 Heliconius species studied differed significantly in wing area and aspect ratio (ANOVA, aspect ratio: F 12, 3502 = 228.4, P < 0.0001, area: F 12, 3502 = 216.4, P < 0.0001; Tukey-adjusted comparisons Fig. S2). We estimated within-species trait repeatability to assess their reliability as species mean estimates for phylogenetic analyses. Wing aspect ratio had higher intra-class repeatability than wing area, with 74% and 48% of the total aspect ratio and size variance explained by differences between species, respectively (aspect ratio: R = 0.74, SE = 0.09, P < 0.0001; size: R = 0.48, SE = 0.1, P < 0.0001). We estimated intraclass repeatability for males and females separately to remove the potential effect of size sexual dimorphism on trait variation, and male size repeatability remained much lower than male wing aspect ratio repeatability (male aspect ratio: R = 0.75, SE = 0.08, P < 0.0001; male size: R = 0.53, SE = 0.1, P < 0.0001). Females had the same wing aspect ratio repeatability as males, whereas wing size repeatability was lower for females probably Mean wing aspect ratio showed no phylogenetic signal (Abouheif's Cmean = 0.15, P = 0.1; Figs. S3 and S5B); in other words, closely related species were not more similar to each other than to distant ones. In contrast, mean wing area showed a strong phylogenetic signal, by which phylogenetically closely related species were more likely to have similar wing areas (Fig. 3, Abouheif's Cmean = 0.33, P = 0.01; Figs. S3 and S6A and B). Wing areas of species in the melpomene clade were on average 14.8% larger than those of species in the erato clade, with H. timareta being 64% larger than H. sara (Fig. 3, H. timareta: mean = 606.6 mm 2 , SE = 3.1; H. sara: mean = 387 mm 2 , SE = 2.9).
Nevertheless, when controlling for sex ratios and allometry on the traits under study, wing aspect ratio and size, the residuals of both traits show a strong phylogenetic signal (Figs. S5 and S6A and C; aspect ratio residuals: Abouheif's Cmean = 0.42, P < 0.001; Fig. S3A and C; size residuals: Abouheif's Cmean = 0.44, P < 0.001). These results support the use of phylogenetic models to study variation in wing aspect ratio and size across species.

PATTERNS ACROSS SPECIES AND ALTITUDES
Species mean altitude had an effect on wing area and aspect ratio (Table 1). Species wings got rounder, that is, lower aspect ratios, with increasing altitudes both when accounting for fixed effects and the phylogeny (Table 1; for full model, see Table S4). These patterns were also evident when examining raw mean wing aspect ratios (Fig. 4A, Gaussian LM: F 1, 9 = 5.37, P < 0.05, R 2 = 0.30), except in the H. telesiphe and H. clysonymus highlands clade, which showed significant phylogenetic autocorrelation (Moran's I index: H. clysonymus 0.53, H. telesiphe 0.49). Species wings got larger with elevation (Table 1; for full model, see Table S4). Without accounting for phylogeny or any fixed effect this is only evident in the erato clade, where high altitude species were bigger than their lowland sister species (Fig. 4B, blue, Gaussian LM: F 1,10 = 17.1, R 2 = 0.80, P = 0.03). However, when assessing individuals from all species together, it becomes clear that larger individuals of both clades tend to be found at higher altitudes (Fig. S8). Both wing size and wing aspect ratio were also significantly correlated with distance from the Equator, and wing aspect ratio was affected by species sex ratio too (Table S4).

PATTERNS WITHIN SPECIES AND ACROSS ALTITUDES
Wings got rounder (lower aspect ratio) with increasing altitude in H. erato and H. melpomene (Fig. 5, H. erato: Gaussian LM: F 6, 1296 = 32.7, P < 0.001, R 2 = 0.13; H. melpomene: Gaussian LM: F 6, 673 = 20.1, P < 0.001, R 2 = 0.14). Individual altitude was the strongest predictor of wing aspect ratio for both species, with sex and wing area being second best in H. erato and H. melpomene, respectively (Table S6 and Fig. S13A and B; Fig. 5). Conversely, the relative importance of explanatory variables of wing area varied for each species (Table S6 and Fig. S13A and B; Fig. 5), and the H. erato model explained less of the overall variation in wing area (Fig. S11, H. erato: Gaussian LM: F 7,1295 = 9.36, P < 0.001, R 2 = 0.04, H. melpomene: Gaussian LM: F 7, 672 = 23.06, P < 0.001, R 2 = 0.18). Wing area in H. erato was correlated with allometric factors interacting with altitude, whereas wing area in H. melpomene was correlated with distance from the Equator (Table S6 and Figs. S10 and S13C and D). Wing area and aspect ratio differed among co-mimicking races of H. erato and H. melpomene, despite inhabiting the same geographic areas (Fig. S12).

Discussion
The fascination for butterfly wing coloration has stimulated many generations of research and Heliconius wing patterns have proven to be excellent study systems for understanding evolution and  H. telesiphe, H. clysonymus, H. erato, H. eleuchia, H. sara, H. doris, H. xanthocles, H. hierax, amH. wallacei, H. numata, H. melpomene, H.   timareta, and H. cydno. Two species, H. telesiphe and H. clysonymus, showed high levels of phylogenetic autocorrelation (Fig. S7) and were thus excluded from the linear model plotted (but not from the main analyses where phylogeny is accounted for).  0.008 * * speciation. Here, we have extended this research by examining wing shape and size variation among more than 3500 individual butterflies, across sexes, clades, and altitudes in 13 species of Heliconius butterflies. We have shown that a large proportion of female biased sexual size dimorphism can be explained by the evolution of larval gregariousness, and that male biased sexual size dimorphism is present only in species that lay eggs individually, regardless of their mating strategy. For the first time in this system, we describe wing morphological variation across environmental clines, with species and populations found at higher altitudes consistently having rounder wings. Here, we demonstrate that Heliconius wing area and aspect ratio are potentially shaped by a plethora of behavioral and environmental selection pressures, in addition to those imposed by Müllerian mimicry.

WING ASPECT RATIO VARIATION
Wing aspect ratio in butterflies and other flying animals determines flight mode and speed (Farney and Fleharty 1969;Buler et al. 2017), and is therefore predicted to vary with life-history requirements across sexes and species. Despite being a simple descriptor of wing shape, aspect ratio has been demonstrated to correlate functionally with gliding efficiency in butterflies by increasing lift-to-drag ratios (Ortega Ancel et al. 2017;Le Roy et al. 2019). Long wings are generally associated with faster gliding flying, whereas round wings with low aspect ratio values favor slow but more maneuverable flight motions (Betts and Wootton 1988;Chai and Srygley 1990;Chazot et al. 2016;Le Roy et al. 2019). For instance, monarch butterfly populations with longer migrations have more elongated wings than resident populations (Satterfield and Davis 2014), and males of Morpho species that dwell in the canopy also have higher aspect ratios to glide faster through open areas (DeVries et al. 2010). In contrast, female Morpho butterflies tend to have rounder wings, and aspect ratio sex differences are stronger in species with colour dimorphism, as varying crypsis may require specific flight behaviors (Chazot et al. 2016).
Heliconius are not notoriously sexually dimorphic especially when compared to other butterflies such as Morpho (Chazot et al. 2016;Jiggins 2016). However, there are important behavioral differences between the sexes. Females are thought to have different flight habits, as they spend much of their time looking for specific host plants for oviposition (Dell'Aglio et al. 2016), or precisely laying eggs on suitable plants, while males tend to patrol open areas searching for receptive females and visit flowers more often (Joron 2005;Jiggins 2016). Thus, it might be predicted that females should have lower aspect ratios, that is, rounder wings, than males (Jones et al. 2013). However, we only found three species with significant, but opposing, sexually dimorphic wing aspect ratios. The wings of males in H. erato were longer than the wings in females, whereas male H. eleuchia and H. wallacei had rounder wings than those of females (Fig. S3). Heliconius wing shape sex differences may require multivariate descriptors of wing morphology and/or analysis of the hindwings, which possess the pheromone-dispersing androconial patch in males (Jones et al. 2013;Mérot et al. 2013Mérot et al. , 2016. In addition, the relatively low collection numbers of female Heliconius could hinder the detection of subtle wing aspect ratio differences across the sexes. Sexual selection has long been known to affect wing color pattern in Heliconius, as it is used for mate recognition and choice (Merrill et al. 2012). More recently, wing aspect ratio has been shown to be part of the mimetic warning signal in Heliconius and their co-mimics (Jones et al. 2013), as it determines flight motion and defines the overall appearance of the butterfly (Srygley 1994(Srygley , 2004a. For instance, wing aspect ratios between two different morphs of H. numata differed consistently across their overlapping ranges, in parallel with their respective and distantly related Melinea co-mimics (Jones et al. 2013). Within-morph wing aspect ratio variation was observed across the altitudinal range of H. timareta in Peru (Mérot et al. 2016), and in the Heliconius postman mimicry ring in Brazil significant across-species wing aspect ratio differences were also found (Rossato et al. 2018a). These studies highlight that while it is clear that color pattern and, to some extent, flight are important for mimicry in Heliconius, wing aspect ratio is also subject to other selection pressures (Mérot et al. 2016;Rossato et al. 2018b).
We found that species inhabiting higher altitudes tend to have rounder wings, after accounting for phylogeny, sample size, and intraspecific variance (Fig. 4A), except in the H. telesiphe-H. clysonymus clade. The latter species may require morphometric analyses of wing tip shape alone, as the overall wing morphology differs significantly from the rest of the Heliconius species here studied (Fig. S7). Interestingly, these patterns were maintained within-species, with high-altitude populations of H. erato and H. melpomene having lower aspect ratios (Fig. 5). Furthermore, altitude was the best predictor of wing aspect ratio in both species (Fig. S13). Rounder wings aid maneuverability and are associated with slower flight in butterflies (Berwaerts et al. 2002;Le Roy et al. 2019) and slower flights are generally associated with a decrease in ambient temperature (Gilchrist et al. 2000). In addition, air pressure, which directly reduces lift forces required to offset body weight during flight (Dillon 2006), decreases approximately 12% across the mean altitudinal range of the species here studied. Thus, the rounder wings in high altitude Heliconius species and populations may aid flying in dense cloud forests, where increased maneuverability could be beneficial, or compensating for lower air pressure at higher altitude.

WING AREA VARIATION
Wing area showed significant sexual dimorphism in more than half of the species studied here, but some species had larger males and others larger females (Fig. 2). In most butterflies, females are overall larger than males, presumably because fecundity gains of increased body size are greater for females (Allen et al. 2011). Larger wings are required to carry larger and heavier bodies, and so Lepidoptera females also tend to have larger wings (Allen et al. 2011;Le Roy et al. 2019). Indeed, in this study, the Heliconius species with larger-winged females were those that lay eggs in large clutches and that have highly gregarious larvae (Fig. 2, Beltrán et al. 2007). A recent study on two species not included here reported wing size dimorphism with larger females in the gregarious H. eratosignis ucayalensis and larger males in the singleegg layer H. demeter joroni (Rosser et al. 2019). Thus, females of these species are likely investing more resources in fecundity than males, which leads to larger body and wing sizes that allow them to carry and lay eggs in clutches throughout adulthood. Larval development time correlates with adult size in H. erato (Rodrigues and Moreira 2002) and growth rates seem to be the same across sexes, at least in the gregarious H. charithonia (Kemp 2019), so we hypothesize that females take longer to develop in gregarious species. Selection for larger females is generally constrained by a trade-off between the benefits of increased fecundity at the adult stage and the higher predation risk at the larval stage associated with longer development times (Allen et al. 2011). This constraint might be alleviated in the unpalatable larvae of Heliconius, as bigger larval and adult size could increase the strength of the warning toxic signal to predators (Jiggins 2016).
An extensive survey identified that only six percent of lepidopteran species exhibit male-biased sexual size dimorphism, and that these patterns were generally explained by male-male competition (i.e., intrasexual selection), in which larger males had a competitive advantage (Stillwell et al. 2010;Allen et al. 2011). In contrast, nearly half of the Heliconius species studied here have male-biased sexual size dimorphism, and all of these lay eggs individually and have solitary larvae (Fig. 2). Male-male competition is high for Heliconius species, as females rarely re-mate despite their very long reproductive life-spans (Merrill et al. 2015). In addition, large reproductive investments in the form of nuptial gifts from males can, in principle, explain male-biased sexual size dimorphisms, as is the case in the polyandrous butterfly Pieris napi whose male spermatophore contains the amount of nitrogen equivalent to 70 eggs (Karlsson 1998;Allen et al. 2011). Male Heliconius spermatophores are not only nutrient-rich, but also loaded with anti-aphrodisiac pheromones that prevent re-mating of fertilized females (Schulz et al. 2008;Merrill et al. 2015). Therefore, it seems likely that in species that lay eggs individually, sexual selection favoring larger males exceeds selection pressures for the large female size needed to carry multiple mature eggs. To our knowledge, Heliconius is the first example of a butterfly genus in which both female-and male-biased size dimorphism are found and can be explained by contrasting reproductive strategies.
We found a strong phylogenetic signal for wing area, with species from the erato clade being on average 12% smaller than those in the melpomene clade (Fig. 3). There are many ecological factors that could explain this pattern, and all could have contributing effects that are hard to disentangle (Fig. 3). First, the erato clade is characterized by facultative pupal-mating (Beltrán et al. 2007;Jiggins 2016), by which males fight for pupae, guard them, and mate with females as they are emerging from the pupal case (Deinert et al. 1994;Jiggins 2016). Smaller males have been shown to outcompete others for a spot on the female pupal case and more successfully inseminate emerging females compared to larger, less agile males (Deinert et al. 1994), which would remove the potential choice of females for larger males. Second, pupal-mating seems to have far-reaching impacts on species lifehistories (Boggs 1981). Species in the melpomene or adult-mating clade are polyandrous, which leads to selection favoring large spermatophores (Boggs 1981) to provide mated females with abundant nutritional resources and defenses that prevent them from re-mating with other males (Cardoso et al. 2009;Cardoso and Silva 2015). This could decrease selection pressure for larger males in the pupal-mating clade, as nuptial gifts need not be so large or nutrient/defense rich, leading to smaller male and female offspring. However, the single origin of pupal-mating in Heliconius (Fig. 2) makes it challenging to infer the impacts of this mating strategy on wing morphology, as the behavior is confounded by phylogeny.
Wing area across species positively correlated with altitude in the erato clade (Fig. 4B), but no clear pattern was found for the melpomene clade species here studied. In contrast, wing area variation within-species (H. erato and H. melpomene) was more correlated with geography (distance to Equator, longitude) and allometry than with altitude (Fig. S10). Nevertheless, high-altitude populations of H. melpomene were slightly bigger than their lowland conspecifics, whereas H. erato did not change (Fig. S13). Two major environmental factors are known to affect insect size across altitudinal clines. One is temperature, such that at lower temperatures, development times are longer and insects grow larger (Chown and Gaston 2010). This perhaps explains cases of Bergmann's rule among ectotherms, where larger species are found in colder climates (Shelomi 2012;Classen et al. 2017). In the geographical range here studied (Fig. 1), we predict temperatures to vary more dramatically along elevation gradients than latitudinal gradients (García-Robledo et al. 2016). We found some evidence for species being bigger with increasing latitudes when accounting for phylogeny and allometry (Table S4), in accordance with Bergmann's rule, but more species at the extremes of the ranges are needed to clarify this (Fig. S7).
Wing beat frequency tends to be lower at low temperatures, so larger wings are required to compensate and gain the extra lift required for flight, as seen in Drosophila robusta (Azevedo et al. 1998;Dillon 2006). A second factor likely to contribute to altitude related differences in wing area is air pressure changes and the correlated lower oxygen availability, which affects flight motion and kinematics as well as many physiological processes. Highaltitude insects can minimize the impacts of lower air pressure by having larger wings, because this lowers the velocity required to induce flight (Dudley 2002).

HERITABILITY
Our study demonstrates that multiple selective forces may be affecting Heliconius wing area and aspect ratio. However, this raises the question of how plastic these traits are in the wild. In Drosophila, the genetic architecture of wing aspect ratio appears to be complex (Gilchrist and Partridge 2001), and is independent of that of wing area (Carreira et al. 2011). Within-species variability of wing area halved when flies were reared in controlled conditions compared to wild populations whereas wing shape variability remained the same, but both traits had a detectable and strong heritable component (Bitner-Mathé and Klaczko 1999). In this study, we found that 74% of the variation in wing aspect ratio could be explained by species identity, in contrast to 48% of the variation in wing area. This high and moderate intraclass repeatability is indicative of heritable traits (Nakagawa and Schielzeth 2013). The fact that closely related species are more likely to have similar wing morphologies, that is, phylogenetic signal, is also indicative of species-level heritability (Queiroz and Ashton 2004).
In insects, wing shape is functionally more constrained than wing size. For example, genetic manipulations of wing shape in Drosophila melanogaster have shown that even subtle changes can have huge biomechanic impacts (Ray et al. 2016), whereas wing/body size differences may impact fecundity more than survival. Here, we find size differences between sexes that can be explained by reproductive strategy, and are likely to be genetically controlled as most sexual dimorphisms are (Allen et al. 2011).The patterns of variation in size across altitudes or latitudes are often not due to phenotypic plasticity alone, as many studies have shown their retention when populations are reared in common-garden conditions (Chown and Gaston 2010). In Monarch butterflies, for example, common-garden reared individuals from wild populations that had different migratory habits showed a strong genetic component for both wing aspect ratio and size (Altizer and Davis 2010).
We have shown that different selection pressures may be shaping the evolution of wing morphology in Heliconius and that the strength of these varies across sexes and environmental clines. Interestingly some of these patterns are maintained at the intraspecific level, with high-altitude populations of H. erato and H. melpomene having rounder wings (Fig. 5), thus potentially adapting locally to the environment in the same way that species of this genus have adapted to altitude over longer evolutionary timescales (Fig. 4). Future work should assess the adaptive significance, plasticity, and heritability of these traits with commongarden rearing and physiological assays in controlled conditions.

Conclusions
Here, we have demonstrated how an understanding of natural and evolutionary history can help to disentangle the putative agents of selection on an adaptive trait. Wing trait differences across sexes, clades, and environments give insight into the selective forces driving phenotypic divergence in Heliconius, beyond the effects of natural selection imposed by Müllerian mimicry. Our study highlights the complexity of selection pressures affecting seemingly simple traits and the need for a thorough understanding of life history differences amongst species.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Table S1. Study species summary data. Table S2. Study species sexual dimorphism. Table S3. Weighted PGLS model selection table for species 15 sexual size dimorphism 16 (SSD), mean wing aspect ratio and mean wing area based on AICc. Table S4. Phylogenetic Generalised Least Squares full 25 model summaries for sexual 26 size dimorphism, wing shape and wing size. Table S5. Model selection based on AIC of within species variation 32 in wing aspect 33 ratio and wing area of H. erato and H. melpomene. Table S6. Full model output table for within-species (39 H. erato and H. melpomene) 40 analyses of wing aspect ratio and wing area. Figure S1. Number of Heliconius individuals in this study collected across 3-year 48 intervals. Figure S2. Wing area (mm2, A) and wing aspect ratio (wing roundness, B) variation 53 across species. Figure S3. Abouheif C-mean distribution plots for six variables. Figure S4. Sexual wing aspect ratio dimorphism across species of the erato cade 69 (A) and the melpomene clade (B). Figure S5. Phylogenetic signal in wing shape. Figure S6. Phylogenetic signal in wing size. Figure S7. Local Moran's I index values for each species for wing area mean (left) 104 and wing aspect ratio mean (right). Figure S8. Wing area variation with altitude across individuals from all species of the 111 erato clade (blue) and the melpomene clade (orange). Figure S9. Species variation in wing area. Plot shows the correlation between 120 distance from the Equator (degrees) and species mean wing area (mm2). Figure S10. Within-species variation in wing area (mm2) across alt.s in H. erato 130 (blue) and H. melpomene (orange), females (left) and males (right). Figure S11. Species variation in raw wing aspect ratio (A) and wing area (B) in H. 140 erato (blue) and H. melpomene (orange). Figure S12. Wing aspect ratio (A) and area (B) variation across mimicry ring wing 149 patterns of the two most abundant species, H. erato (blue) and H. melpomene 150 (orange). Figure S13. Relative importance of model predictors of within species variation wing 157 aspect ratio (A, B) and wing area (C, B) in H. erato (A, C) and H. melpomene (B, D).