Effects of natural wing damage on flight performance in Morpho butterflies: what can it tell us about wing shape evolution?

Understanding the evolution of wing morphology requires estimating the impact of wing shape variation on fitness-related behaviours. In butterflies, the capacity to fly enabled by wing morphology is crucial throughout adult life during many key behaviours involved in survival, such as resource foraging (Hall and Willmott, 2000) or escaping from predators (Barber et al., 2015; Chai and Srygley, 1990), and in reproduction, such as male–male contest (Berwaerts et al., 2006; Wickman, 1992) or courtship (Scott, 1974). Wing shape directly affects various aspects of flight performance, ranging from the energy budget (Ancel et al., 2017) to the aerodynamic forces produced during wingbeats (Ellington, 1984; Muijres et al., 2017). Investigating the consequences of wing shape variation on these different components of flight performance may shed light on the forces driving wing shape evolution within and across species (Arnold, 1983; Norberg and Rayner, 1987; Le Roy et al., 2019). For example, selection acting on wing shape has been evidenced by comparing migrating and non-migrating populations of monarch butterflies: migrating individuals exhibit a more elongated wing shape, probably reducing flight cost and hence benefiting their long-distance migration (Altizer and Davis, 2010; Dockx, 2007). However, the consequences of wing shape variation on flight performance are poorly documented, preventing the precise identification of the selective pressures acting on wing shape evolution (Chazot et al., 2016; Johansson et al., 2009; Outomuro and Johansson, 2015).

To investigate the effect of wing shape variation on flight performance, most previous studies have involved experimental manipulation of wing shape. For example, experimental wing clipping in butterflies has highlighted the importance of the hindwings for linear and turning acceleration in Lepidoptera (Jantzen and Eisner, 2008). Artificially modified wings have been used to investigate how insects compensate for changes in wing morphology by altering their behaviour. The loss of wing surface has been shown to induce an increase of wingbeat frequency in bees (Hedenström et al., 2001; Vance and Roberts, 2014), flies (Muijres et al., 2017) and butterflies (Kingsolver, 1999). In hawkmoths, asymmetrical wing loss causes the insect to flap its damaged wing with a larger amplitude, correcting for the unequal force production (Fernández et al., 2012). In flies, compensation for asymmetrical wing loss involves a body roll towards the damaged wing and changes in wing motion (Muijres et al., 2017). Although behavioural adjustments in response to wing damage may have evolved in some insects, a significant fitness cost of wing damage has nevertheless been documented in some studies. Field studies on bumblebees have shown a reduced foraging efficiency in damaged individuals (Higginson and Barnard, 2004) as well as a lower life expectancy (Cartar, 1992). In dragonflies, reduction of hindwing area significantly decreases capture success (Combes et al., 2010). In butterflies, however, a mark–recapture study found no effect of wing surface reduction on survival (Kingsolver, 1999).

In the studies mentioned above, wing alterations were typically generated by clipping the wings, and usually consisted of a symmetric versus asymmetric reduction of the wing surface (Fernández et al., 2012, 2017; Haas and Cartar, 2008; Vance and Roberts, 2014), or a gradual reduction in the wing surface (Muijres et al., 2017; Vance and Roberts, 2014). But these artificial alterations may not reflect the natural wing damage experienced by insects in the wild. The spatial distribution and extent of damage on the wings of wild insects have rarely been quantified in natural populations (but see Higginson and Barnard, 2004). In the wild, a frequent source of wing damage is collision with obstacles (Foster and Cartar, 2011). Accidental collisions may rip the wing in various degrees and directions, although this is somewhat constrained by the vein architecture (Mountcastle and Combes, 2014; Wootton, 1992). Collisions may also occur during agonistic interactions with conspecifics (Alcock, 1996; Carvalho et al., 2016). Failed predator attacks can also cause significant damage in different wing locations (Carpenter, 1942; Edmunds, 1974; Robbins, 1981; Shapiro, 1974). As a result, flying insects have to cope with a wide diversity of damage during their life, varying in both extent and location. While some damage may have too harsh a consequence on flight to even be observable in the wild, most frequently observed damage is probably less deleterious, or more easily compensable. Studying wing damage in wild-caught individuals that probably has a limited impact on survival capacity should allow characterization of wing shape variation with a diversity of effects on fitness, from neutral to severely deleterious. This situation is therefore relevant to the investigation of the consequences on flight of a biologically realistic range of damage, as a step towards inferring the selective pressures acting on wing shape.

In this study, we used natural and artificially accentuated damage on butterfly wings to test how modification of different wings and different wing areas affects flight performance. We studied wild-caught individuals of two species of large butterflies, Morpho helenor and Morpho achilles (Nymphalidae, Satyrinae), exhibiting a large variation in wing damage. We first quantified their flight performance using three-dimensional videography. We then precisely determined the frequency and spatial location of damage over the four wings using geometric morphometrics, allowing us to estimate and compare the effects of natural and artificially accentuated damage on flight performance. Although the initial cause of wing damage in the wild is not known, our study captures part of the range of naturally occurring wing shape alteration. The reported effects on flight performance are thus expected to reflect typical challenges that butterflies face in the wild throughout their adult life. By generating more severe wing damage in our sample, we aimed at contrasting natural – presumably non-lethal – damage, with damage of higher magnitude that would more strongly threaten survival in the wild.

Field sampling was performed in July and August 2016, in the middle basin of the Río Huallaga (San Martín Department, Peru), near the city of Tarapoto, along the Rio Shilcayo (06°27′07″S, 76°20′47″W; ca. 450 m a.s.l.) and the village of San Antonio de Cumbasa (06°24′24″S, 76°24′25″W; ca. 470 m a.s.l.). We sampled a total of 63 Morpho, including 32 individuals from the species Morpho helenor (Cramer 1776) and 31 from the species Morpho achilles (Linnaeus 1758). These two species are nearly identical phenotypically (Blandin, 1988) and are thought to have similar flight behaviours (DeVries et al., 2010). Only three females were caught, largely because of their cryptic flight behaviour, contrasting with the extensive patrolling displayed by males. Half of the captured specimens were undamaged, while the rest had their wings at least slightly damaged. In order to increase the variation of wing damage in our sample, 25 of the captured specimens were stored in a mesh-cage (4 m×2 m×1.8 m) for 3 days so as to generate collision-induced wing damage. The final sample was composed of 31 undamaged individuals and 32 individuals for which the extent of damage ranged from a small amount to half of the wing surface missing. Out of the 32 damaged individuals, 13 came from the wild and 19 from the cage.

Butterflies were filmed in a large outdoor insectary (8 m×4 m×2.5 m) built close to the sampling sites, in a sheltered spot where no wind was detectable. Each specimen was released from a shaded side of the cage, and generally flew towards the sunniest part of the cage. Two video cameras (GoPro Hero4 Black set at 60 images per second) mounted on a tripod at fixed height were used to record the films. In order to capture most of the flight path, the camera zoom lens was set on wide angle (focal length: 14 mm), thereby allowing coverage of the entire volume of the cage when combining the two camera views. A successful sequence was defined as a flight path moving through the entire field of view of the two cameras. Multiple trials were performed until a minimum of three successful sequences were obtained for each individual. We recorded a total of 227 successful sequences with a mean duration of 1.4±0.89 s, ranging from 0.6 to 6.8 s. After being filmed, specimens were placed in a −20°C freezer, ensuring subsequent morphological measurements exactly matched the state in which butterflies had flown.

Sequences of the same flight obtained from the two cameras were first synchronized with respect to a reference frame. The frame distortion (fisheye effect) due to wide-angle settings was corrected without limiting the wide view angle, using the DWarp Argus package (Jackson et al., 2016) implemented in Python 2.7. To recover exact distances from films, cameras were calibrated with the direct linear transformation (DLT) technique (Hartley and Sturm, 1995; http://www.kwon3d.com/theory/dlt/dlt.html) by digitizing an object of known length (here a wand) moved throughout the experimental cage. The wand tracking and DLT coefficient computation were performed using DLTdv5 (Hedrick, 2008) and easyWand (Theriault et al., 2014) Matlab program, respectively. To facilitate the tracking of the filmed butterflies, a background subtraction algorithm (KaewTraKulPong and Bowden, 2002) was applied to each video via Python 2.7. Trajectory points were digitized for each frame at the centroid of the butterfly. Wingbeats were quantified by manually digitizing a point on the butterfly within frames containing the highest upstroke positions and those containing the lowest downstroke positions transcribing the spatial and temporal position of each wing stroke along the flight trajectory (Fig. 1). Based on the temporal wing stroke position, gliding flight phases along the trajectory were defined as at least 10 consecutives frames (representing 0.16 s) without any wing stroke. The other parts of the trajectory were considered to be flapping phases (Fig. 1). By reducing wing surface, wing damage can limit the gliding capacity of the observed butterflies: this can be behaviourally compensated by a reduction in the length of the gliding phases and/or an increase in wing beat frequency during the flapping phases. To distinguish these two types of behavioural compensation, we computed both the total wing beat frequency through time and the wing beat frequency within the flapping phases.

Flight trajectories were then smoothed using a low-pass Butterworth filter (order 4, critical frequency 0.5) (Butterworth, 1930), removing the high-frequency movements, i.e. the steep gaps – artefactual movements within trajectories, stemming from tracking noise (see Fig. S1). Using a custom-written R script, we calculated the following parameters for each flight trajectory: flight velocity, sinuosity (computed as the ratio of the actual distance covered over the distance between start and end position), wingbeat frequency, gliding proportion and flight height. These non-redundant parameters were chosen for their relevance in describing flight behaviour. They were computed as the mean value over the flight trajectory. We also extracted the maximal duration of gliding and flapping phases, and the smallest turning angle of each trajectory as a measure of manoeuvrability (Rayner, 1988). We summarized the flight performance of each individual by retaining the mean and maximum values out of the three flight sequences analysed per individual.

The forewings (FWs) and hindwings (HWs) were photographed in dorsal view using a Nikon D90 camera in controlled light conditions. The spatial variation of damage within and between wings, but also between naturally and artificially damaged individuals, was visually assessed by generating a heat map (Fig. 2; see also Figs S2 and S3). This was done by superimposing standardized images of the wings from all specimens using EBImage (R package; Pau et al., 2010) so as to count the missing area at the pixel scale. The heat map was built from a matrix summing the occurrences of missing pixels, and plotted with autoimage (R package; French, 2017).

As our goal was to capture the different distribution of damage throughout the wings, we considered damage variation as a wing shape alteration, and quantified it using a landmark-based geometric morphometric method (Adams et al., 2004; Bookstein, 1997). This method has proven to be well suited for studying variation in butterfly wing shape and venation (e.g. Breuker et al., 2010; Chazot et al., 2016; Zhong et al., 2016), although it has never been used to quantify variation in wing damage. Because vein architecture differed between specimens depending on wing damage, we did not use it to place fixed landmarks. Instead, we used 300 semi-landmarks equidistantly spaced along the (more or less damaged) wing outline. Semi-landmarks can be used when identifiable landmarks are unavailable. To remove the variation along the outline due to a lack of homology, the semi-landmarks are allowed to slide along the local tangent to the curve in an iterative process (Gunz and Mitteroecker, 2013). Once slid, they can be treated as regular landmarks in the analyses. We also placed one fixed landmark at the wing base, fixing the overall landmark configuration with respect to this homologous position available for all specimens. The procedure was applied to both the left- and right-reflected FWs and HWs. All landmarks were digitized using TpsDig2 (Rohlf, 2015).

In order to obtain variables describing shape alteration between specimens, we then performed a generalized procrustes analysis (Rohlf and Slice, 1990) on the landmark configurations of each of the four wings using geomorph (R package; Adams and Otárola-Castillo, 2013). This procedure extracts the shape information from landmark positions by removing the extraneous variations; namely, the position, size and orientation (Adams et al., 2004). Variation in the newly obtained procrustes coordinates then only reflects shape variation. Because our goal was to relate wing shape alteration to that measured in flight parameters, it was necessary to analyse simultaneously both pairs of wings. Indeed, the observed flight results from the combined use of the four wings, and is therefore potentially affected by their combined shape alteration. After each wing was superimposed separately, we conducted a principal component analysis (PCA) on each sets of procrustes coordinates. We then combined the principal components (PCs) accounting for 90% of cumulated variance of each of these four PCAs. Finally, we performed a new PCA on this global dataset to obtain PCs that combine the shape information of the four wings per individual (Fig. S4). We also applied this procedure independently on the FWs and HWs to focus on the effect of their respective shape alteration on flight separately.

Wing area was computed from the previously digitized landmarks using Momocs (R package; Bonhomme et al., 2014) independently on each of the four wings. Wing area was preferred over other measures of size (e.g. centroid size): it is indeed directly relevant for aerodynamics and was found to be less affected by extreme irregularities in outlines. For specimens with only the left (or right) wing damaged, the extent of the damage was calculated by subtracting the area of the damaged wing from the area of its intact counterpart. When specimens had both wings damaged, we used the mean wing area of the intact specimens. Damage extent was expressed separately for the FWs and HWs to test their respective effect on flight performance. We also computed the asymmetry resulting from wing damage, expressed as the absolute value of the difference in the extent of damage between the right and left sides of the specimen.

All statistical analyses were conducted in R (http://www.R-project.org/). We first tested the effect of species and sex on flight parameters, while considering only undamaged specimens, using multivariate analysis of variance (MANOVA). Variation in flight performance was examined by performing a PCA on the set of flight parameters (Fig. S5 and Table S1). PCs are linear combinations of the original data that maximize the variation between individuals. The first PCs, accounting for the major trends in flight variation, were used as flight data in the subsequent multivariate analyses. The co-variation between wing damage and flight was first assessed by Escoufier's RV coefficient (Escoufier, 1973; Klingenberg, 2009). We then performed a two-blocks partial-least squares (2B-PLS) regression between flight and shape datasets. The 2B-PLS analysis specifically focuses on the co-variation between two sets of multivariate data, by constructing pairs of linear combinations of variables within each dataset (here, the shape and flight PCs) that maximally co-vary across datasets (Rohlf and Corti, 2000). Then, in order to identify the precise effect of area loss on flight behaviour, we tested the effect of FW and HW area loss and their asymmetry on each flight parameter, using multiple regression analysis. Species, sex and cause of damage (i.e. cage or nature) were included as factors to control for their respective effects on flight.

To test for the effects of the distribution of damage on flight performance, we performed multiple regressions of the shape PCs on each flight parameter. Finally, because the PCs are composite variables combining information on different wings, direct visualization of the associated morphologies is not straightforward. To visualize damage variation explained by the different PCs, we separated individuals into three groups of equal size depending on their location along the PCA axis (first, second and last third), and used the heat maps to depict damage variation within each group (see below).

The extent of damage within our sample ranged from 0 to 39.49% of missing wing area at the individual level, with up to 45.96% and 51.46% for the FWs and HWs, respectively. Wing asymmetry was strongly correlated with extent of damage (r=0.72; P<0.001), highlighting that wing damage rarely occurred symmetrically. The proportion of damage between FW and HW was not correlated (r=0.15; P=0.31) but both wings were on average similarly affected (mean extent of damage 13.37±1.90% on FWs versus 10.68±1.75% on HWs; P=0.37, W=1222). Individuals kept in a mesh-cage showed greater wing damage relative to those damaged in nature. This difference was mostly due to greater damage on the FWs in individuals kept in a mesh-cage (P<0.001, W=0.76 versus P=0.05, W=616 when comparing FW and HW damage, respectively). It should be noted, however, that wing damage produced in the cage may have added to natural damage, for which we have no record. Types of damage were nevertheless different between individuals from the cage and those from nature. Collisions occurring in the cage mostly affected the leading edge of the FWs. In contrast, the most frequent natural damage was located along the marginal and posterior zone of the HW (Fig. 2).

As similar results were found when using either the mean or the maximal values of flight parameters, only the maximal values are considered below. Amongst all flight parameters computed, only gliding proportion and longest gliding phase were strongly correlated (r>0.50): we thus excluded the longest gliding phase from the analysis. The MANOVA performed on undamaged individuals revealed no difference among species on flight parameters (Wilks' λ=0.854; P=0.52). Although our sample included only two undamaged females, an effect of sex was detected on flight speed: female flight was significantly slower than male flight (P<0.05, W=54).

Co-variation between wing shape and flight parameters was significant (RV coefficient=0.20; P=0.01). The first PLS vector (53.87% of covariation explained) carried a variation in flight parameters (PLS 1 block 2) opposing higher flight speed and gliding proportion to longer flapping phase duration and higher wingbeat frequency (Table 1). This was associated with the shape component (PLS 1 block 1) describing an accumulation of wing damage mostly located on the FWs (Fig. 3). This co-variation was clearly driven by severely damaged individuals, as both intact and slightly damaged individuals displayed similar variation in flight parameters. Another pattern of co-variation was detected on the second PLS axis (21.33% of co-variation explained), where most of the flight variation was driven by the flight height, and associated with variation in wing damage located on the HWs.

When running multiple regression analysis, we found that the effect of wing area loss on flight parameters was not affected by sex (P=0.91) or by the cause of damage (i.e. cage versus nature; P=0.50). Although no effect of species was detected on flight parameters among undamaged individuals (M. helenor: n=14; M. achilles: n=17), a significant effect of species was detected on flight speed (P<0.05) and gliding proportion (P<0.05) when testing for the effect of wing damage on flight. Controlling for the species effect, we found that FW area loss had a significant negative effect on flight speed (R²=0.26, F_60,0.35=10.69, P<0.001) and gliding proportion (R²=0.19, F_60,24.09=7.10, P<0.001). Specifically, M. achilles and M. helenor both flew more slowly when their FWs were damaged, although for the same degree of damage, M. helenor reduced its normal flight speed by 19% (0.43 m s⁻¹) while M. achilles reduced it by only 13% (0.32 m s⁻¹). Similarly, the extent of gliding flight decreased as FW damage increased for both species, but M. achilles reduced its gliding proportion by 59% while M. helenor only reduced it by 32% for a same degree of damage (Fig. 4). A slight increase in wingbeat frequency was associated with FW area loss (R²=0.11, F_60,1.2=3.67, P<0.05), but such an effect was not detected when focusing only on flapping phases (R²=0.05, F_60,2.1=1.63, P=0.50), consistent with a transition from flap–gliding flight to continuous flapping flight associated with FW damage. FW area loss was indeed positively associated with maximal flapping duration (R²=0.23, F_61,0.57=18.67, P<0.001). HW area loss did not contribute to the variation in these flight parameters, but had a significant negative effect on flight height (R²=0.08, F_61,0.26=5.65, P<0.05). No effect of wing asymmetry was detected, probably because it was strongly correlated with wing area loss.

The different types of damage distribution captured by the shape variation of the PCs had contrasting effects on flight parameters. In particular, the reduction in flight speed and the extension of flapping phase duration were explained by FW shape alteration described on PC2 (R²=0.18, F_61,0.37=13.61, P<0.001 and R²=0.22, F_61,0.58=17.71, P<0.001, respectively). These two flight parameters were negatively correlated (r=−0.41; P<0.001): most damaged individuals both had a reduced flight speed and used almost only flapping flight. In contrast to the FW shape variation carried on PC2, variation on PC1 showed no relationship with flight variation, neither did any PC carrying HW shape variation. The distribution of damage throughout the wings shown by the heat maps revealed that damage located at the FW tips (shape variation on PC2) was associated with reduced flight speed. Damage occurring along the FW margin, in contrast, showed no effect on variation in flight speed (Fig. 5).

Our results show that highly damaged M. helenor and M. achilles display a reduction in the typical flap–gliding flight observed in intact or less damaged individuals, progressively replaced by strict flapping flight. As continuous flapping flight is more energetically costly than gliding flight (Dudley, 2002), the reduced gliding ability in very damaged individuals may result in increased flight cost, ultimately impacting survival and/or reproductive abilities.

In previous studies investigating wing damage in insects, an increase in wingbeat frequency in damaged individuals has frequently been reported (Fernández et al., 2012; Hedenström et al., 2001; Kingsolver, 1999; Muijres et al., 2017; Vance and Roberts, 2014), often associated with increased metabolic costs. This kinematic adjustment allows maintenance of sufficient lift in spite of a reduced wing area (Altshuler et al., 2005; Dickinson et al., 1998). Higher wingbeat frequency following wing loss was measured during hovering flight in the butterfly Pontia occidentalis (Kingsolver, 1999) and the moth Manduca sexta (Fernández et al., 2012). Damaged Morpho butterflies showed a slight increase in wingbeat frequency, at least during the forward flights studied here. This increased frequency mostly stemmed from the smaller number of gliding phases observed in damaged individuals, because the wingbeat frequency during flapping phases did not significantly increase in damaged individuals as compared with intact ones. In our experiment, Morpho butterflies thus did not modulate wingbeat frequency in response to wing damage, but mostly limited gliding phases, switching from a flap–gliding flight to more continuous flapping, with similar wingbeat frequency. It has been shown that reduced lift can be balanced by an increase in either wingbeat frequency or stroke amplitude (Altshuler et al., 2005). Morpho may thus also compensate for wing damage by adjusting wing stroke amplitude rather than wingbeat frequency. A more precise comparison of kinematic parameters (such as wing stroke amplitude) of intact and damaged individuals would be needed to test this hypothesis. Changes in wingbeat frequency in response to wing damage may also depend on the type of flight muscle involved. In asynchronous flyers (such as flies or bees) a single nervous impulse triggers multiple contractions, allowing a high wingbeat frequency to be attained (Dudley, 1991). In contrast, synchronous flyers (such as butterflies or dragonflies) have a one to one correspondence between nervous impulses and muscle contraction (Pringle, 1981). Physiological differences among insect species might impact their ability to adjust wingbeat frequency. Wing clipping in synchronous flyers indeed produces little effect on wingbeat frequency (Roeder, 1951). The increase in wingbeat frequency associated with clipped wings was found to be approximately 2 Hz in moths (Fernández et al., 2017), while an increase of 19 Hz and 23 Hz was measured in bees and flies, respectively (Muijres et al., 2017; Vance and Roberts, 2014). Whether these changes in frequency reflect a behavioural adjustment or a passive mechanical response is unknown. The capacity to increase wingbeat frequency may nonetheless vary among insects, leading to different constraints on the evolution of wing shape and toughness. Morpho butterflies showed a mean wingbeat frequency of 5.9±1.7 Hz, ranking them among the lowest frequencies found in insects (Sotavalta, 1947). Such a low wingbeat frequency is expected given that their flight is composed of frequent alternations between flapping and gliding phases.

Although rarely tested, how the distribution of missing area over the wings impacts flight behaviour would pinpoint the contribution of different wing parts to flight performance. Our results clearly suggest an unequal impact of FW and HW damage on Morpho flight: damage on the FWs more strongly affected flight performance than that on the HWs. The loss of FW area significantly reduced flight speed and the proportion of gliding flight. Experimental wing manipulations have shown that complete removal of the FWs makes butterflies flightless, while they can still fly without the HWs (Jantzen and Eisner, 2008). Interestingly, we found that damage specifically located at the tip of the FWs – altering the shape of the leading edge – most strongly impaired the flight speed and gliding ability of Morpho butterflies, probably because the leading edge of the wing is strongly involved in lift generation: during flight, the incoming air flow separates precisely as it crosses the leading edge, producing a vortex (termed the leading edge vortex) that creates a suction force resulting in lift enhancement (Ellington et al., 1996; Sane, 2003). Butterflies with a deteriorated leading edge face a substantial reduction in lift generation.

While the most anteriorly located wing damage had severe consequences on flight, the loss of other wing parts seemed to be relatively harmless: expansion of damage located along the FW margin had no detectable effect on flight. Damage on the HWs showed only a limited effect on flight parameters: the clearest impact was on flight height, which was reduced in damaged individuals, indicating a possible role for the HWs in upward flight. Such an effect needs to be tested in a more controlled experiment. In particular, it would be interesting to investigate the role of the HWs for upward escape from predators, especially during take-off. This could have important consequences for HW shape evolution. Jantzen and Eisner (2008) showed that HW removal in butterflies was associated with a significant reduction in linear and turning acceleration, limiting performance in zig-zag, erratic flight, and therefore putatively decreasing the capacity to escape flying predators. Here, we did not find such an effect, possibly because the distribution of damage over the four wings varied greatly between individuals, limiting our statistical power.

Morpho helenor and M. achilles are sister species showing extreme phenotypic similarity and occupying the same microhabitat (Blandin, 1988; Chazot et al., 2016). Nevertheless, our results show that behavioural changes triggered by wing damage differed between these two Morpho species. While wing area loss resulted in a decrease in flight speed and gliding proportion for both species, flight speed decreased significantly more in M. helenor than in M. achilles. In contrast, the gliding ability of M. achilles was more impaired by wing damage relative to that of M. helenor. However, no difference in any flight parameters was found between the two species when considering only intact individuals, and similar damage proportions were observed in the two species. This behavioural difference may stem from subtle differences in the location of lost areas among damaged M. helenor and M. achilles, differently impacting their flight. Nevertheless, by constraining butterflies to fly with severely damaged wings, we may also have revealed differences in flight behaviour or morphology between species (such as muscles mass, position and power) that would otherwise have remained undetectable (i.e. in less challenging conditions). Challenging conditions eliciting extreme performance have indeed been shown to reveal the consequences of morphological or physiological variations more readily than favourable conditions (Losos et al., 2002; Wainwright and Reilly, 1994). This apparent difference between M. achilles and M. helenor in their behavioural ability to compensate for wing damage suggests the need to consider the interactions between wing shape and other behavioural and morphological traits when investigating the evolution of wing shape across butterfly species.

Our results highlight that impaired flight performance (and possible behavioural compensation) depends not only on the extent of the damage but also on its location over the wings. The crucial role of the leading edge in flight, for instance, might generate a strong selection on its toughness. The evolution of a close proximity and even fusion of several veins in the anterior part of insect wings might also stem from such selection on wing toughness (Dudley, 2002; Rees, 1975). A deteriorated leading edge was indeed rarely observed in our sample of naturally damaged individuals, and was only found in captive individuals experiencing frequent collisions with the rough cage walls. This extreme damage contrasts with that observed in wild-caught individuals, such as scattered tearing along the wing margin, with limited impact on flight performance. The damage observed in wild-caught individuals probably illustrates a range of wing shape alterations with limited impact on survival. Such harmless consequences of wing margin damage may explain the evolution of eyespots along the wing margin in some butterfly species (e.g. in Bicyclus anynana), which have been shown to deflect predator attacks away from vital body parts (Lyytinen et al., 2003; Stevens, 2005). Further field studies quantifying natural wing damage in butterflies should assess the frequency of scattered margins in the wild, particularly in species displaying peripheral eyespots. Experimental manipulations of wing shape are still required to rigorously test the effect of quantitatively similar but spatially different wing loss, to identify the selection regimes affecting the various parts of the wing and thus altering wing shape evolution.

Because alterations of FW shape have a much more severe impact on flight performance, FW shape may be under stronger stabilizing selection than HW shape. Strauss (1990) found that wing shape variation in Heliconiinae and Ithomiinae butterflies increased from the anterior part of the FW to the posterior part of the HW. Such an increase in shape variability and conversely a decreasing density of veins along the chordwise wing section may reflect an antero-posterior decrease in aerodynamic constraints. Lepidoptera HWs frequently show large shape variation between species, such as the presence of scalloped edges or expanded tails (Robbins, 1981; Rubin et al., 2018), contrasting with the generally subtler variation in FW shape. HW tails that closely resemble a butterfly's head in some lycaenid species are thought to deflect predator attacks (López-Palafox and Cordero, 2017; Robbins, 1981). Similarly, tails in moths were recently shown to have a deflecting role against bat attacks (Barber et al., 2015; Rubin et al., 2018). The evolution of marked HW extensions in response to predation suggests that a large shape variation may occur on butterfly HWs with limited effects on flight performance. Aerodynamic constraints on HWs may thus be slighter than those acting on FWs, limiting the deleterious impact of the evolution of HW variations on flight performance.

Altogether, by studying wing shape variation induced by natural damage, our work suggests that contrasted selective regimes may act on different wing parts of Morpho butterflies, highlighting wing areas under stabilizing and relaxed selection. Ascertaining the variation in aerodynamic constraints within and between insect wings may thus provide important insight on the evolution of wing morphology. Further experimental studies generating a large diversity of wing damage should enable a better understanding of these constraints.

The authors would like to thank the Peruvian authorities, and in particular SERFOR (the Servicio Nacional Forestal y de Fauna Silvestre) for providing the necessary research permits (permit: 002-2015-SERFOR-DGGSPFF). We thank Ronald Mori-Pezo, Corentin Clerc, and César, Jay and Bruno Ramirez for help with the sampling of butterflies. We are grateful to Tyson Hedrick for his advice on the three-dimensional videography method. We thank two anonymous reviewers for helpful comments on the first draft of the manuscript.