Genomics of altitude-associated wing shape in two tropical butterflies

Understanding how organisms adapt to their local environment is central to evolution. With new whole-genome sequencing technologies and the explosion of data, deciphering the genomic basis of complex traits that are ecologically relevant is becoming increasingly feasible. Here we study the genomic basis of wing shape in two Neotropical butterflies that inhabit large geographical ranges. Heliconius butterflies at high elevations have been shown to generally have rounder wings than those in the lowlands. We reared over 1100 butterflies from 71 broods of H. erato and H. melpomene in common-garden conditions and show that wing aspect ratio, i.e. elongatedness, is highly heritable in both species and elevation-associated wing shape differences are maintained. Genome-wide associations with a published dataset of 666 whole genomes from across a hybrid zone, uncovered a highly polygenic basis to wing shape variation in the wild. We identify several genes that have roles in wing morphogenesis or wing shape variation in Drosophila flies, making them promising candidates for future studies. There was little evidence for molecular parallelism in the two species, with only one shared candidate gene, nor for a role of the four known colour pattern loci, except for optix in H. erato. Thus, we present the first insights into the heritability and genomic basis of within-species wing shape in two Heliconius species, adding to a growing body of evidence that polygenic adaptation may underlie many ecologically relevant traits.


Introduction
Climate change is forcing organisms to 'move, adapt, or die'. With temperatures rising and land-use changing in the lowlands, shifting to higher elevations might be the only way to flee extinction for many taxa (I.-C. Chen, Hill, Ohlemüller, Roy, & Thomas, 2011). However, the environment changes drastically along mountains, with diverse sets of challenges expected to drive local adaptation (Halbritter, Billeter, Edwards, & Alexander, 2015). Thus, identifying the genomic mechanisms that allow organisms to inhabit wide ranges is key to understanding which taxa are most likely to adapt locally when forced to colonise new, highelevation, environments. With novel sequencing technologies and the explosion of genomic data, we now have the tools to decipher the genetic basis of ecologically relevant traits in the wild.
Insect flight has many essential functions, including dispersal, courtship, and escaping predators (Dudley, 2002). As such, it is under strong selection to be optimised to suit local environments. Air pressure decreases with elevation, which in turn reduces lift forces required for taking flight, as well as oxygen available for respiration (Hodkinson, 2005).
Subtle variation in wing morphology can have big impacts on flight mode and performance (Berwaerts, Van Dyck, & Aerts, 2002). For example, butterflies in tropical rainforests have been found to have rounder wings when they inhabit the understory, both within closely related taxa (Chazot et al., 2014) and across many species (Mena, Kozak, Cárdenas, & Checa, 2020). Elongated wings reduce wing-tip vortices, resulting in more efficient and faster flight, whereas short and wide wings are associated with higher manoeuvrability and lift (Le Roy, Debat, & Llaurens, 2019a). Monarch butterflies have evolved two wing phenotypes; migratory populations have longer wings for long-distance gliding than those that remain year-round in Caribbean islands (Altizer & Davis, 2010). Thus, wing shape represents a trait likely undergoing strong selection across elevations and conferring local adaptation in insects.
Insect wing shape is phylogenetically conserved in many taxa (Houle, Mezey, Galpern, & Carter, 2003;Montejo-Kovacevich et al., 2019;Thomas, Walsh, Wolf, McPheron, & Marden, 2000), but also highly evolvable in the laboratory (Houle et al., 2003), suggesting that wing shape can be heritable. Wing shape tends to have higher heritability than wing size in flies and wasps (Gilchrist & Partridge, 2001;Klingenberg, Debat, & Roff, 2010;Xia, Pannebakker, Groenen, Zwaan, & Bijma, 2020), as size is more prone to life-history trade-offs and condition-dependency. Temperature changes have been shown to trigger plastic responses in a wide range of traits (Deutsch et al., 2008;Merilä and Hendry, 2014; but see Fraimout et al., 2018), including heat tolerance across altitudes in Heliconius (Montejo-Kovacevich et al., 2020). Thus, it is important to ascertain the heritability of morphological traits when studying their genomic basis, as plasticity could affect trait values and reduce the power to detect genomic associations. Experimental common-garden designs can help overcome the challenges of studying phenotypic clines in the wild, by providing the same environmental conditions to genotypes of different populations, which allows for the estimation of heritability (de Villemereuil, Gaggiotti, Mouterde, & Till-Bottraud, 2016) Wing morphology has been mostly studied in Drosophila, with developmental pathways identified for D. melanogaster (Carreira, Soto, Mensch, & Fanara, 2011;Diaz de la Loza & Thompson, 2017) and many of the candidate genes functionally tested (Pitchers et al., 2019;Ray, Nakata, Henningsson, & Bomphrey, 2016). In other insects, wing dimorphisms have been particularly well-studied. These tend to have simple genetic architectures controlling wingless, or short winged, morphs (B. Li et al., 2020;McCulloch et al., 2019), or are environmentally induced with hormonal regulation, such as in migratory planthoppers (Xu & Zhang, 2017;Zhang, Brisson, & Xu, 2019). However, studies describing the genetic basis of quantitative wing shape variation in organisms other than Drosophila are lacking.
Significant advances have been made in understanding the genetic basis of local adaptation in the wild. There are two general strategies to identify loci potentially under selection: (i) forward genetics, where known phenotypic traits are associated to genotypes (via genomewide association studies or quantitative trait loci with laboratory crosses), and (ii) reverse genetics, where variation in allele frequencies in natural populations is studied to detect signatures of selection across the genome, without any prior knowledge of the phenotypes involved (Fuentes-Pardo & Ruzzante, 2017;Pardo-Diaz, Salazar, & Jiggins, 2015). On the one hand, population genomics, i.e. reverse genetics, excels at detecting genomic regions under selection across environments, but the lack of phenotypic associations can make findings hard to interpret. On the other hand, GWA studies, i.e. forward genetics, require quantifiable traits on a large number of individuals, thus making it difficult to study populations in the wild.
A good solution is to study steep clines, where the environment changes continuously over a small space while gene flow is high, and combine both forwards and reverse genetics approaches (Cornetti & Tschirren, 2020;Tigano & Friesen, 2016). Sequencing individuals along such clines allows for sufficient phenotypic variance to measure and associate to genotypes (forward genetics), while maintaining low genetic structure, and additionally test which genomic regions might be undergoing selection by comparing the extremes of the clines (reverse genetics).
The study of aposematic wing pattern colouration in Heliconius butterflies is a prime example of this approach. Decades of careful lab crosses described the inheritance of the colour pattern elements and identified the genes controlling them (Martin et al., 2012;Nadeau et al., 2016;Reed et al., 2011;Westerman et al., 2018). In the last decade, whole-genome sequencing in elevationally structured colour-morph hybrid zones has confirmed the role of these loci, which repeatedly diverge across highland and lowland colour pattern races (reverse genetics, Meier et al., 2020;Nadeau et al., 2014). Genome-wide association studies demonstrated that these diverging parts of the genome were associated with specific colour pattern elements. These loci have now been functionally tested in Heliconius and other taxa, making important contributions to our understanding of evolution and speciation. Thus, whole-genome sequencing populations across steep environmental clines, while studying phenotype heritability with controlled rearing, is a good approach to disentangle the genomic underpinnings of ecologically relevant traits.
Here we study the genetic basis of wing shape variation in two widespread species of Heliconius butterflies across an elevational cline in the Ecuadorian Andes. Heliconius inhabiting high altitudes have recently been found to have rounder wings than lowland butterflies, a pattern seen both across species and within species along elevational clines . To estimate the heritability of this potentially adaptive trait, we common-garden reared 71 broods of H. erato lativitta and H. melpomene malleti from across the cline, yielding 1141 offspring (Fig. 1). We then used forward (GWAS) and reverse (Fst) genetic approaches with whole-genome data of 666 individuals to identify regions associated with quantitative variation in wing aspect ratio and determined which regions diverged between extremes of the cline. This genomic dataset was obtained from a study that developed a new low-cost linked-read sequencing technology, 'haplotagging', to examine colour pattern clines in an altitudinally structured hybrid zone (Fig. 1C , Meier et al., 2020). Here we present the first study to examine the heritability and genomic basis of wing shape of two butterfly species in the wild.  . B) Common-garden rearing protocol. C) Topographic surface of transect across elevations used for whole-genome sequencing. Both species cooccur and have three main colour pattern morphs along this cline: two distinct colour pattern morphs (H. e notabilis and H. m. plesseni, referred to as "black", and H. e lativitta and H. m. malletti, referred to as "rayed") and within-species hybrids displaying admixed phenotypes (green circles), the most common hybrid phenotype is depicted.

STUDY SYSTEM AND WILD BUTTERFLY COLLECTION
H. erato and H. melpomene are the two most widespread Heliconius species, and have Müllerian aposematic mimicry to advertise their toxicity to predators (Chris D. Jiggins, 2016).
They can be found continuously coexisting across elevational clines ranging from sea level up to 1600 m along the Andean mountains. We sampled females of H. erato lativitta and H. melpomene malleti across the eastern slope of the Ecuadorian Andes for common-garden rearing (Fig. 1B, orange triangles). For the genomic analyses we used a large dataset from a nearby hybrid zone

WING MEASUREMENTS
Wing morphology was analysed with an automated pipeline in the public software Fiji (Schindelin et al., 2012), using custom scripts from Montejo-Kovacevich et al., 2019. These automatically crop, extract the right or left forewing and perform particle analysis on the wing ( Fig. 1B). Wing area is estimated for the whole wing in mm 2 . For wing aspect ratio, the particle analysis function obtains the best fitting ellipse of the same area as the wing, extracts the major and minor axis' lengths, from which aspect ratio is estimated (major axis/minor axis). We only include forewings in this study, as they determine flight speed and mode, whereas hindwings act as an extended surface to support flight and gliding (Le Roy, Debat, & Llaurens, 2019b;Wootton, 2002). Furthermore, hindwings tend to be more damaged in Heliconius due to in-flight predation and fragile structure.

COMMON-GARDEN REARING
Fertilised females of H. erato lativitta and H. melpomene malleti were caught in the wild with hand nets at elevations ranging from 380 m up to 1600 m (Table S1). Females from all altitudes were simultaneously kept in separate 2x1x3 m cages of purpose-built insectaries at the Universidad Regional Amazónica Ikiam (

WHOLE-GENOME DATASET
We used 666 whole-genomes of H. erato (n=479) and H. melpomene (n=187) from a recent study , sequenced with 'haplotagging'. This linked-read sequencing technique retains long-range information via barcoding of DNA molecules before sequencing, which permits megabase-size haplotypes to be reconstructed computationally . More details on analysis and phasing of molecules can be found in the Supplementary Materials of Meier et al., 2020, and summarised

STATISTICAL ANALYSES
All non-genomic analyses were run in R V2.13 (R Development Core Team 2011) and graphics were generated with the package ggplot2 (Ginestet, 2011). Packages are specified below and all R scripts are publicly available (Zenodo: TBC). Sequence data from  is deposited at the NCBI Short Read Archive (PRJNA670070).

Effects of altitude on wing shape
To test the effects of maternal altitude on wing aspect ratio of common-garden reared offspring, we fitted a linear mixed model that included as fixed effects wing area, sex, development time (days from larva hatching until pupating), and altitude ("high" if the mother was collected above 600m in elevation) with lme4 model fits (Bates, Mächler, Bolker, & Walker, 2015). All continuous fixed effects were standardized to a mean of zero and unit variance to improve model convergence (Zuur, Ieno, Walker, Saveliev, & Smith, 2009). We included family ID as a random effect (intercept) to account for relatedness among offspring in H. melpomene. In H. erato, we nested family ID within experiment location as an additional random effect, as 10 of the families were reared in common-garden insectary conditions (2018), whereas the rest were reared in laboratory conditions (2019). We performed backward selection of random and fixed effects, in that order, with the package lmerTest (Kuznetsova, Brockhoff, & Christensen, 2017), with likelihood ratio tests and a significance level of α = 0.1 (functions ranova() and drop1() Kuznetsova et al., 2017). When comparing models with different fixed effects we fitted Maximum Likelihood (REML = FALSE), otherwise Restricted Maximum Likelihood models were fitted. Model residuals were checked for homoscedasticity and normality. With the coefficient of determination (R 2 ), we estimated the proportion of variance explained by the fixed factors alone (marginal R 2 , R 2 LMM(m)) and by both, the fixed effects and the random factors (conditional R 2 , R 2 LMM(c)), implemented with the

Heritability estimates
To test the heritability of wing shape, we assessed wing aspect ratio variation across individuals from families reared in common-garden conditions. We used all broods with at least three offspring that could be phenotyped. First, we first used an ANOVA approach, with family identity as a factor explaining the variation in aspect ratio. We then estimated narrowsense heritability (h 2 ) with two approaches (more details can be found in Note S2). Firstly, we estimated intra-class correlation coefficient (ICC) or repeatability with a linear mixed model approach; family ID was set as the grouping factor, with a Gaussian distribution and 1000 parametric bootstraps to quantify uncertainty, implemented with the function rptGaussian() in the rptR package (Stoffel, Nakagawa, & Schielzeth, 2017). Secondly, we estimated narrowsense heritability (h 2 ) from the slope of mother and mid-offspring wing aspect ratio regressions for those families where the mother's wings were intact (31/48 broods in H. erato and 10/23 in H. melpomene). Father's phenotypes were unknown but their effect should be random, hence unbiased, with respect to the mother's phenotypes.

Genome-wide association mapping of wild wing aspect ratio
All population genomics analyses were performed in ANGSD version 0.933, which uses genotype likelihoods as input to account for genotype uncertainty and a Bayesian framework well-suited for large low-coverage sequencing datasets (Korneliussen, Albrechtsen, & Nielsen, 2014). To account for population structure across the cline, we first calculated admixture proportions with NGSadmix. We used genotype likelihoods computed by STITCH as input and ran NGSadmix on a random subset of 10% of the total sites with a minor allele frequency of at least 0.05, and specified two (k=2) or three (k=3) clusters for H. erato and H. melpomene, respectively.
To identify the genomic regions potentially controlling quantitative wing shape variation in these two Heliconius species, we performed genome-wide association mapping (GWAS, doAsso ANGSD). VCF files from Meier et al. (2020) containing genotype likelihoods were used as input (-vcf-pl). We performed the GWAS with a generalized linear framework and a dosage model (-doAsso 6), which calculates the expected genotype from the input and implements a normal linear regression with the dosages as the predictor variable. Aspect ratio was used as the continuous response variable (-yQuant). Three covariates were incorporated (-cov) to control for sex, wing area, and population structure with the admixture proportions obtained from NGSadmix, as described above. Likelihood ratio test (LRT) statistics are calculated per site (following a chi square distribution with one degree of freedom) under an additive model with logistic regression (-model 1) (Skotte, Korneliussen, & Albrechtsen, 2012).
In GWAS with fewer unlinked SNPs, those with strong associations can sometimes be found in isolation, i.e. without flanking SNPs showing association, and assumed to be in linkage with the causal SNP. However, with large SNP datasets that are not pruned for LD, many are expected to show strong associations when located near the causal site or at random due to noise (Zhou et al., 2020). To partly account for this, we obtained median and minima SNP association p-values for sliding windows of 50 SNP, with a step size of 10 SNP, with the R package 'WindowScanR' (Tavares, 2016(Tavares, /2020. When visualising our results, we first use median p-value per window, plotted as -log10(p-value), to regionally smooth the results so that spurious associations get dampened by their flanking high p-value SNPs, while regions with many SNPs with strong associations will be easily identifiable, as their median will remain high. This is analogous to recently developed methods for medical genetics that use  (Hinrichs, 2006).

Identifying regions diverging between highland and lowland populations
We computed genetic differentiation (Fst) between the highland and lowland subspecies (excluding the mid-elevation hybrids), to identify regions diverging across altitudes and potentially overlapping with wing shape associated regions. With the genotype likelihoods estimated by STITCH, we calculated the site frequency spectrum (SFS) for each population (dosaf, ANGSD). Then, we obtained folded 2D-SFS for both populations combined to use as a prior for the joint allele frequency probabilities with the function realsfs. Fst was calculated per site using the Weir-Cockerham correction (realsfs fst index, ANGSD), and 5kb window averages with 1kb steps were obtained for plotting (following Meier et al. 2020).

Candidate gene annotation
Regions associated with wing shape that had at least one 50 SNP window with a median pvalue below the threshold of p<0.01 (-log10(p)median>2), and ten windows with lower or equal slope=0.79, R 2 =0.6**). Due to the strong sexual wing shape dimorphism in H. erato, motherto-male offspring regressions had a higher intercept (SI Fig. S2A). There was no correlation between mother's and offspring's wing size in H. erato but a strong correlation in H. melpomene (Fig. S3, H.

Figure 2.
Mother and mid-offspring regression for wing aspect ratio (A,C) and F1 offspring wing aspect ratio with respect to maternal origin across elevations (B,D). Each black point represents mean wing aspect ratio per family and its size is proportional to number of offspring per family. Grey shading around the regression corresponds to 95% confidence intervals of the regression. Stars represent significance levels of two sample t-tests between highland (>600 m) and lowland (<600 m) family means for each species (+<0.075, *< 0.05, **<0.01, ***<0.001).

Subspecies and hybrid phenotypes
Aspect ratio varied slightly across the elevational cline sampled for whole-genome sequencing of both species (Fig. 1 C). Highland subspecies, H. e. notabilis and H. m.
plesseni (black, Fig. 3) were on average rounder than hybrids (green, Fig. 3), and in H. erato they were also rounder than lowland subspecies H. e. lativitta (orange, Fig. 3 A). Hybrids did not differ in wing aspect ratio compared to their corresponding lowland subspecies, in other words, hybrids were phenotypically more lowland-like (Fig. 3A, Fig. S7, S8). More importantly, the aspect ratio of hybrid individuals encompassed most of the trait variance of the pure subspecies (Fig. S6), increasing our power to detect genomic associations.

Association mapping of aspect ratio variation
Genome-wide association mapping revealed a highly polygenic basis for wing shape variation in both species (Fig. 4) . We obtained association statistics for 11.3M SNPs ( Fig. 4, Fig. S9). Despite a highly polygenic basis, there is evidence that certain genomic regions were strongly associated (Q-Q plots Fig. S10, permutations Fig. S11). We further investigated 28 regions that had 10 adjacent outlier windows supported by permutations and had a median p-value<0.01 (-log10(p)>2, Fig. 4).

Candidate genes
We identified genes in the 28 regions of interest, within or neighbouring the strongest SNP hits, yielding 23 and 22 candidate genes for H. erato and H. melpomene, respectively. In 12 out of the 28 regions, SNPs within the region of strongest association were found within genes that could be annotated (Table 2, 3). The remaining outliers were intergenic (n=16) and potentially associated with regulatory variants, and were, on average, 35.8kb from the nearest (n=6) or second nearest gene (n=10, Of the known loci that control colour pattern differences across this cline (Fig. S12, , only one was significantly associated with wing shape variation. The transcription factor controlling presence/absence of red, optix, was strongly associated with wing shape in H. erato, whereas in H. melpomene only a few windows were associated, thus not meeting our criteria for 'region of interest' (Fig. 5B, Fig. S11). The outlier region on chromosome 13 was found in close proximity across species, with both SNPs with the lowest p-value mapping nearby the rugose gene, which affects neuromuscular junction development, synaptic architecture, brain morphology, and associative learning (Fig. 5A,  Zhao et al., 2013). The epidermal growth factor receptor (Egfr) regulates rugose (Shamloula et al., 2002) and has been associated with wing shape variation in wild D.

FST GENOME SCANS
There is little background genomic differentiation between highland and lowland populations of both species (mean FST in H. erato: 0.0261 and in H. melpomene: 0.0189, Fig. S12, Meier et al., 2020). Of the 28 regions identified as potentially associated with wing shape variation, ten were found to be in regions of elevated genomic differentiation between highland and lowland populations (FST green lines, Fig. 4), from low levels of differentiation (e.g. H. erato 1(b), Fstmax=0.08) to high differentiation (e.g. H. erato 13, Fstmax=0.44). The strongest four FST peaks in this cline are associated with colour patterning (Fig. S9, Meier et al., 2020;Nadeau et al., 2014), but only optix (chr. 18) in H. erato was also strongly associated with wing shape (Fig. 5B).

Discussion
Here we combine the power of hybrid zones across steep environmental clines, common garden rearing, and whole-genome sequencing to study the genomic basis of a potentially adaptive trait in the wild. We found that wing aspect ratio is highly correlated between mothers and their offspring in two butterfly species, highly repeatable across commongarden reared offspring families, and correlated with the altitude at which the mother was collected (Fig. 2). With a large dataset comprising 666 whole-genomes sequenced with haplotagging  and association mapping, we uncover a highly polygenic basis to wing shape, and identify potential candidate genes in regions with many SNPs showing associations (Fig. 4). Furthermore, with a population genetics approach, we find that many of these regions are also diverging between highland and lowland populations, potentially being selected for local adaptation to highland environments.

WING SHAPE IS HERITABLE
The amount of wing shape variation explained by family across common-garden reared offspring was high for both species (H. erato: 21% and H. melpomene: 39%), especially when compared to the 74% of variance explained by species identity in a previous comparative study ( Montejo-Kovacevich et al., 2019). The resemblance in wing aspect ratio between mothers and their offspring is indicative of a highly heritable trait (Fig. 2 A), although we cannot rule out maternal effects partly driving this pattern. Rearing offspring in commongarden conditions strongly reduces the effects of shared mother-offspring environmental variables, but cannot account for, for example, variation in resources the mothers provide to their eggs. We found, however, that mother wing size did not correlate with offspring wing sizes in H. erato, whereas it did in H. melpomene. This highlights that wing shape might be less affected by maternal effects compared to other more condition-dependent traits, such as size. Furthermore, strong sexual dimorphism in wing shape present in H. erato, was maintained in common-garden reared individuals and both sexes had similar correlations with mother phenotypes, implying a strong genetic component to wing shape variation (Allen, Zwaan, & Brakefield, 2011).
From a local adaptation perspective, we might predict wing shape to be highly heritable.
Insects can behaviourally compensate for damaged or abnormal wings through changes in flight and body kinematics (Fernández, Driver, & Hedrick, 2017). Yet, this might incur a fitness cost, as many behaviours, such as courtship and predator escape, are dependent on efficient flight (Le Roy et al., 2019a). Generally, cases of wing shape plasticity are rarer than size plasticity, especially if the two traits are allometrically decoupled, allowing for subtle changes to be selected if advantageous (Carreira et al., 2011;Gilchrist & Partridge, 2001;Strauss, 1990). In Heliconius, wings have been found to be rounder at higher elevations, both across and within species that inhabit large ranges (Montejo-Kovacevich et al., 2019). In our study, wing shape differences observed in the wild in H. erato and H. melpomene were maintained in common-garden reared broods, with individuals from highland mothers having, on average, rounder wings (Fig. 2B, Fig. S4). Together, this supports the hypothesis that subtle changes in wing aspect ratio are highly heritable, and may be involved in local adaptation to altitude.

OTHER INSECTS
We found 5/28 regions mapping to genes involved in the biological process of 'wing disc development', which is remarkable given that there are only 419 genes under this category out of 64,000 genes on Flybase (Thurmond et al., 2019), plus one involved in wing vein formation (Lunde et al., 2003). In H. erato, the most promising candidate gene was the suppressor of deltex, su(dx) (Fig. 4 A, Table 2), an E3 ubiquitin-protein ligase of the Notch signalling pathway (Jennings, Blankley, Baron, Golovanov, & Avis, 2007). su(dx) knockouts in D. melanogaster result in rounder wings via reduction of longitudinal wing venation (Mazaleyrat et al., 2003), wing margin reduction (Wilkin et al., 2004) or via interactions with other proteins (Djiane et al., 2011). Interestingly, this region has moderate levels of genetic differentiation across high and low elevation populations (Fstmax=0.24, Fig. 4A, 1(a)), which could be an indication of altitude-associated selection on this candidate gene.
The knirps-related protein (knrl) is involved in second wing vein development (Table 3, Lunde et al., 2003), and Tap42 (Fig. 4B) triggers apoptosis in the developing wing discs (Wang et al., 2012). Interestingly, a recent study found that knockdowns of fat and knirps affected wing shape in D. melanogaster, which regulate two of our candidate genes (Table 3, Pitchers et al., 2019). Here we show a highly polygenic basis to aspect ratio variation in two Heliconius species, which contrasts to the simple genetic architecture of wing dimorphisms in other insects (B. Li et al., 2020). The additive effects of many genes could facilitate altitudeassociated adaptations, by refining wing shapes to suit the local environment (Wellenreuther & Hansson, 2016). However, aspect ratio may correlate to other wing shape descriptors, thus future studies could focus on multivariate shape variation in these species.

NOVEL CANDIDATE GENES ASSOCIATED WITH WING SHAPE VARIATION
We found 5/28 regions mapping to genes involved in the biological process of 'open tracheal system development', which only has 277 genes on Flybase (Thurmond et al., 2019). Four of these were involved in septate junction assembly function, a category with 83 genes (Thurmond et al., 2019). varicose, an H. erato candidate gene, is essential to septate junction formation, and Drosophila mutants lead to defective tracheal systems and shrivelled wings (Moyer & Jacobs, 2008;V. M. Wu et al., 2007). The veins determine the architecture of the wing, while circulatory and tracheal systems affect the elasticity and physiological functioning of the wing during flight (Pass, 2018). A recent study has shown that butterfly circulatory and tracheal systems are essential for avoiding overheating of the wings by prompting behavioural responses (Tsai et al., 2020). Heliconius are generally slow gliding flyers, which display their warning colour patterns to predators, so their wings could rapidly overheat from prolonged exposure to the sun. The potential role of tracheae and wing venation in Heliconius flight and thermoregulation remains to be explored.

CONVERGENCE AND COLOUR PATTERN ASSOCIATION
Despite phenotypic convergence towards rounder wings at high altitude, we found little evidence for molecular parallelism underlying wing shape variation between H. erato and H. melpomene. One of the 28 regions identified as potentially involved with this trait was found in the regulatory region of the gene rugose in both species (Fig. 5 A). Mutants of rugose in Drosophila lead to the 'rough eye phenotype', similarly to another candidate in H. melpomene, punch, which results in cone cell loss and defective vision (Shamloula et al., 2002). Rugose mutants also exhibit aberrant associative odour learning, changes in brain morphology, and increased synapses in the larval neuromuscular junction (Volders et al., 2012). A recent study has found Egfr, which regulates rugose, to affect wing shape in Drosophila (Pitchers et al., 2019), making it an interesting candidate for future functional studies. The lack of abundant molecular parallelism in this phenotype is in stark contrast to colour pattern loci in Heliconius, which have been repeatedly co-opted or shared via adaptive introgression to create the diversity of mimetic colour morphs we observe across the Neotropics (Fig. S12, Jiggins, 2016;Nadeau et al., 2014).
We found a strong association with wing shape variation in H. erato at the optix locus, which controls most of the red colour patterning (Bainbridge et al., 2020;Lewis et al., 2019;Meier et al., 2020;Van Belleghem et al., 2017). In Heliconius, wing shape has been traditionally studied in the context of mimicry, as similar wing shapes could aid locomotor mimicry (Mérot, Le Poul, Théry, & Joron, 2016;Srygley, 1994Srygley, , 1999. For example, morphs of H. numata that mimic the distantly related genus Melinea tend to converge in wing shapes where they coexist (Jones et al., 2013). Here, wing shape in H. erato was more variable between highland and lowland subspecies than in H. melpomene (Fig. 3A), which could explain its association with a colour pattern locus (Fig. 5B). It is important to note that, although H. erato and H.
melpomene mimic each other across their range and both have rounder wings in the highlands, the differences in aspect ratio between the two species are larger than those across elevations (Fig. 3D). It therefore seems unlikely that this subtle wing shape variation is primarily caused by mimicry. Furthermore, other species not part of this mimicry ring have been shown to also have rounder wings at high altitudes .
Thus, while some mimicry-related wing shape variation may be controlled by optix, many other loci are likely involved in shaping wings to suit the local environment and life-history of each species.

GENETICS BASIS OF AN ECOLOGICALLY RELEVANT TRAIT
A common criticism of population genomics approaches, reverse genetics, that aim to link genotypes and environments is that they often lack phenotypes. Traits directly measured from the wild might be a result of phenotypic plasticity, and are thus rarely used to infer local adaptation. Common-garden rearing can bridge the gap between phenotypes, genotypes, and environment, by providing measurements of heritability and repeatability of a trait across families whose genetic material comes from different extremes of an environmental cline (de Villemereuil et al., 2016). On the other hand, using highly differentiated populations can result in spurious phenotypic associations and makes identifying divergent outliers challenging. Thus, GWAS in the wild should use randomly mating populations with little population structure, while ensuring there is enough phenotypic variation to detect genetic associations with the trait of interest. Hybrid zones, where closely related subspecies or morphs come into contact along an environmental cline, can provide such ideal conditions to carry out GWAS in the wild.
Here we have demonstrated the value of combining these approaches to disentangle the genomic basis of an ecological relevant trait in the wild. We found that wing aspect ratio is highly heritable in two widespread species of Heliconius butterflies, and that altitude explains part of the variation in this trait. We have identified several regions potentially shaping wings in H. erato and H. melpomene, including five candidate genes involved in wing morphogenesis and several identified to be affecting wing shape in recent Drosophila studies (Pitchers et al., 2019). We found evidence of molecular parallelism between species at the gene rugose, and a strong association of H. erato wing shape with a known colour pattern locus, optix. Our study adds to a growing body of evidence showing that most quantitative traits conferring local adaptation are highly polygenic (Barghi, Hermisson, & Schlötterer, 2020). Spatial environmental heterogeneity and gene flow are thought to maintain high levels of standing genetic variation (Tigano & Friesen, 2016). This can favour polygenic adaptation, so that incomplete sweeps of many redundant loci can shift traits towards an optimum (Yeaman, 2015). A slow-moving optimum, such range-expansions towards the highlands, should favour polygenic adaptation via small-effect loci, whereas selection for a distant optimum, such a switch in colour pattern mimicry in Heliconius, should favour large-effect loci (Barghi et al., 2020). New whole-genome sequencing technologies could foster the study of polygenic adaptation to the environment and shape our understanding of the mode and tempo of evolution in the wild.