Hybrid speciation driven by multilocus introgression of ecological traits.

Hybridization allows adaptations to be shared among lineages and may trigger the evolution of new species1,2. However, convincing examples of homoploid hybrid speciation remain rare because it is challenging to demonstrate that hybridization was crucial in generating reproductive isolation3. Here we combine population genomic analysis with quantitative trait locus mapping of species-specific traits to examine a case of hybrid speciation in Heliconius butterflies. We show that Heliconius elevatus is a hybrid species that is sympatric with both parents and has persisted as an independently evolving lineage for at least 180,000 years. This is despite pervasive and ongoing gene flow with one parent, Heliconius pardalinus, which homogenizes 99% of their genomes. The remaining 1% introgressed from the other parent, Heliconius melpomene, and is scattered widely across the H. elevatus genome in islands of divergence from H. pardalinus. These islands contain multiple traits that are under disruptive selection, including colour pattern, wing shape, host plant preference, sex pheromones and mate choice. Collectively, these traits place H. elevatus on its own adaptive peak and permit coexistence with both parents. Our results show that speciation was driven by introgression of ecological traits, and that speciation with gene flow is possible with a multilocus genetic architecture.

Subtle Introgression Footprints at the End of the Speciation Continuum in a Clade of Heliconius Butterflies.

Quantifying gene flow between lineages at different stages of the speciation continuum is central to understanding speciation. Heliconius butterflies have undergone an adaptive radiation in wing color patterns driven partly by natural selection for local mimicry. Color patterns are also known to be used as assortative mating cues. Therefore, wing pattern divergence is considered to play a role in speciation. A corollary is that mimicry between closely related species may be associated with hybridization and interfere with reproductive isolation. Here, we take a multifaceted approach to explore speciation history, species boundaries, and traits involved in species differentiation between the two closely related species, Heliconius hecale and Heliconius ismenius. We focus on geographic regions where the two species mimic each other and contrast this with geographic regions where they do not mimic each other. To examine population history and patterns of gene flow, we tested and compared a four-population model accounting for linked selection. This model suggests that the two species have remained isolated for a large part of their history, yet with a small amount of gene exchange. Accordingly, signatures of genomic introgression were small except at a major wing pattern allele and chemosensing genes and stronger in the mimetic populations compared with nonmimetic populations. Behavioral assays confirm that visual confusion exists but that short-range cues determine strong sexual isolation. Tests for chemical differentiation between species identified major differences in putative pheromones which likely mediate mate choice and the maintenance of species differences.

The interplay of local adaptation and gene flow may lead to the formation of supergenes.

Supergenes are genetic architectures resulting in the segregation of alternative combinations of alleles underlying complex phenotypes. The co-segregation of alleles at linked loci is often facilitated by polymorphic chromosomal rearrangements suppressing recombination locally. Supergenes are involved in many complex polymorphisms, including sexual, colour or behavioural polymorphisms in numerous plants, fungi, mammals, fish, and insects. Despite a long history of empirical and theoretical research, the formation of supergenes remains poorly understood. Here, using a two-island population genetic model, we explore how gene flow and the evolution of overdominant chromosomal inversions may jointly lead to the formation of supergenes. We show that the evolution of inversions in differentiated populations, both under disruptive selection, leads to an increase in frequency of poorly adapted, immigrant haplotypes. Indeed, rare allelic combinations, such as immigrant haplotypes, are more frequently reshuffled by recombination than common allelic combinations, and therefore benefit from the recombination suppression generated by inversions. When an inversion capturing a locally adapted haplotype spreads but is associated with a fitness cost hampering its fixation (e.g. a recessive mutation load), the maintenance of a non-inverted haplotype in the population is enhanced; under certain conditions, the immigrant haplotype persists alongside the inverted local haplotype, while the standard local haplotype disappears. This establishes a stable, local polymorphism with two non-recombining haplotypes encoding alternative adaptive strategies, that is, a supergene. These results bring new light to the importance of local adaptation, overdominance, and gene flow in the formation of supergenes and inversion polymorphisms in general.

The evolution and ecology of multiple antipredator defences.

Prey seldom rely on a single type of antipredator defence, often using multiple defences to avoid predation. In many cases, selection in different contexts may favour the evolution of multiple defences in a prey. However, a prey may use multiple defences to protect itself during a single predator encounter. Such "defence portfolios" that defend prey against a single instance of predation are distributed across and within successive stages of the predation sequence (encounter, detection, identification, approach (attack), subjugation and consumption). We contend that at present, our understanding of defence portfolio evolution is incomplete, and seen from the fragmentary perspective of specific sensory systems (e.g., visual) or specific types of defences (especially aposematism). In this review, we aim to build a comprehensive framework for conceptualizing the evolution of multiple prey defences, beginning with hypotheses for the evolution of multiple defences in general, and defence portfolios in particular. We then examine idealized models of resource trade-offs and functional interactions between traits, along with evidence supporting them. We find that defence portfolios are constrained by resource allocation to other aspects of life history, as well as functional incompatibilities between different defences. We also find that selection is likely to favour combinations of defences that have synergistic effects on predator behaviour and prey survival. Next, we examine specific aspects of prey ecology, genetics and development, and predator cognition that modify the predictions of current hypotheses or introduce competing hypotheses. We outline schema for gathering data on the distribution of prey defences across species and geography, determining how multiple defences are produced, and testing the proximate mechanisms by which multiple prey defences impact predator behaviour. Adopting these approaches will strengthen our understanding of multiple defensive strategies.

Selective sweeps on novel and introgressed variation shape mimicry loci in a butterfly adaptive radiation.

Natural selection leaves distinct signatures in the genome that can reveal the targets and history of adaptive evolution. By analysing high-coverage genome sequence data from 4 major colour pattern loci sampled from nearly 600 individuals in 53 populations, we show pervasive selection on wing patterns in the Heliconius adaptive radiation. The strongest signatures correspond to loci with the greatest phenotypic effects, consistent with visual selection by predators, and are found in colour patterns with geographically restricted distributions. These recent sweeps are similar between co-mimics and indicate colour pattern turn-over events despite strong stabilising selection. Using simulations, we compare sweep signatures expected under classic hard sweeps with those resulting from adaptive introgression, an important aspect of mimicry evolution in Heliconius butterflies. Simulated recipient populations show a distinct 'volcano' pattern with peaks of increased genetic diversity around the selected target, characteristic of sweeps of introgressed variation and consistent with diversity patterns found in some populations. Our genomic data reveal a surprisingly dynamic history of colour pattern selection and co-evolution in this adaptive radiation.

Dominance mechanisms in supergene alleles controlling butterfly wing pattern variation: insights from gene expression in Heliconius numata.

Loci under balancing selection, where multiple alleles are maintained, offer a relevant opportunity to investigate the role of natural selection in shaping genetic dominance: the high frequency of heterozygotes at these loci has been shown to enable the evolution of dominance among alleles. In the butterfly Heliconius numata, mimetic wing color variations are controlled by an inversion polymorphism of a circa 2 Mb genomic region (supergene P), with strong dominance between sympatric alleles. To test how differences in dominance observed on wing patterns correlate with variations in expression levels throughout the supergene region, we sequenced the complete transcriptome of heterozygotes at the prepupal stage and compared it to corresponding homozygotes. By defining dominance based on non-overlapping ranges of transcript expression between genotypes, we found contrasting patterns of dominance between the supergene and the rest of the genome; the patterns of transcript expression in the heterozygotes were more similar to the expression observed in the dominant homozygotes in the supergene region. Dominance also differed among the three subinversions of the supergene, suggesting possible epistatic interactions among their gene contents underlying dominance evolution. We found the expression pattern of the melanization gene cortex located in the P-region to predict wing pattern phenotype in the heterozygote. We also identify new candidate genes that are potentially involved in mimetic color pattern variations highlighting the relevance of transcriptomic analyses in heterozygotes to pinpoint candidate genes in non-recombining regions.

The gene cortex controls mimicry and crypsis in butterflies and moths.

The wing patterns of butterflies and moths (Lepidoptera) are diverse and striking examples of evolutionary diversification by natural selection. Lepidopteran wing colour patterns are a key innovation, consisting of arrays of coloured scales. We still lack a general understanding of how these patterns are controlled and whether this control shows any commonality across the 160,000 moth and 17,000 butterfly species. Here, we use fine-scale mapping with population genomics and gene expression analyses to identify a gene, cortex, that regulates pattern switches in multiple species across the mimetic radiation in Heliconius butterflies. cortex belongs to a fast-evolving subfamily of the otherwise highly conserved fizzy family of cell-cycle regulators, suggesting that it probably regulates pigmentation patterning by regulating scale cell development. In parallel with findings in the peppered moth (Biston betularia), our results suggest that this mechanism is common within Lepidoptera and that cortex has become a major target for natural selection acting on colour and pattern variation in this group of insects.

When Do Opposites Attract? A Model Uncovering the Evolution of Disassortative Mating.

AbstractDisassortative mating is a rare form of mate preference that promotes the persistence of polymorphism. While the evolution of assortative mating and its consequences for trait variation and speciation have been extensively studied, the conditions enabling the evolution of disassortative mating are still poorly understood. Mate preferences increase the risk of missing mating opportunities, a cost that can be compensated by a greater fitness of offspring. Heterozygote advantage should therefore promote the evolution of disassortative mating, which maximizes the number of heterozygous offspring. From the analysis of a two-locus diploid model with one locus controlling the mating cue under viability selection and the other locus coding for the level of disassortative preference, we show that heterozygote advantage and negative frequency-dependent viability selection acting at the cue locus promote the evolution of disassortative preferences. We predict conditions of evolution of disassortative mating coherent with selection regimes acting on traits observed in the wild. We also show that disassortative mating generates sexual selection, which disadvantages heterozygotes at the cue locus, limiting the evolution of disassortative preferences. Altogether, our results partially explain why this behavior is rare in natural populations.

Inversion breakpoints and the evolution of supergenes.

The coexistence of discrete morphs that differ in multiple traits is common within natural populations of many taxa. Such morphs are often associated with chromosomal inversions, presumably because the recombination suppressing effects of inversions help maintain alternate adaptive combinations of alleles across the multiple loci affecting these traits. However, inversions can also harbour selected mutations at their breakpoints, leading to their rise in frequency in addition to (or independent from) their role in recombination suppression. In this review, we first describe the different ways that breakpoints can create mutations. We then critically examine the evidence for the breakpoint-mutation and recombination suppression hypotheses for explaining the existence of discrete morphs associated with chromosomal inversions. We find that the evidence that inversions are favoured due to recombination suppression is often indirect. The evidence that breakpoints harbour mutations that are adaptive is also largely indirect, with the characterization of inversion breakpoints at the sequence level being incomplete in most systems. Direct tests of the role of suppressed recombination and breakpoint mutations in inversion evolution are thus needed. Finally, we emphasize how the two hypotheses of recombination suppression and breakpoint mutation can act in conjunction, with implications for understanding the emergence of supergenes and their evolutionary dynamics. We conclude by discussing how breakpoint characterization could improve our understanding of complex, discrete phenotypic forms in nature.

How do predators generalize warning signals in simple and complex prey communities? Insights from a videogame.

The persistence of distinct warning signals within and between sympatric mimetic communities is a puzzling evolutionary question because selection favours convergence of colour patterns among toxic species. Such convergence is partly shaped by predators' reaction to similar but not identical stimulus (i.e. generalization behaviour), and generalization by predators is likely to be shaped by the diversity of local prey. However, studying generalization behaviour is generally limited to simple variations of prey colour patterns. Here, we used a computer game played by humans as surrogate predators to investigate generalization behaviours in simple (4 morphs) and complex (10 morphs) communities of unprofitable (associated with a penalty) and profitable butterflies. Colour patterns used in the game are observed in the natural populations of unprofitable butterfly species such as Heliconius numata. Analyses of 449 game participants' behaviours show that players avoided unprofitable prey more readily in simple than in complex communities. However, generalization was observed only in players that faced complex communities, enhancing the protection of profitable prey that looked similar to at least one unprofitable morph. Additionally, similarity among unprofitable prey also reduced attack rates only in complex communities. These results are consistent with previous studies using avian predators but artificial colour patterns and suggest that mimicry is more likely to evolve in complex communities where increases in similarity are more likely to be advantageous.

Hybridization and transgressive exploration of colour pattern and wing morphology in Heliconius butterflies.

Hybridization can generate novel phenotypes distinct from those of parental lineages, a phenomenon known as transgressive trait variation. Transgressive phenotypes might negatively or positively affect hybrid fitness, and increase available variation. Closely related species of Heliconius butterflies regularly produce hybrids in nature, and hybridization is thought to play a role in the diversification of novel wing colour patterns despite strong stabilizing selection due to interspecific mimicry. Here, we studied wing phenotypes in first- and second-generation hybrids produced by controlled crosses between either two co-mimetic species of Heliconius or between two nonmimetic species. We quantified wing size, shape and colour pattern variation and asked whether hybrids displayed transgressive wing phenotypes. Discrete traits underlain by major-effect loci, such as the presence or absence of colour patches, generate novel phenotypes. For quantitative traits, such as wing shape or subtle colour pattern characters, hybrids only exceed the parental range in specific dimensions of the morphological space. Overall, our study addresses some of the challenges in defining and measuring phenotypic transgression for multivariate traits and our data suggest that the extent to which transgressive trait variation in hybrids contributes to phenotypic diversity depends on the complexity and the genetic architecture of the traits.

Polymorphism at a mimicry supergene maintained by opposing frequency-dependent selection pressures.

Explaining the maintenance of adaptive diversity within populations is a long-standing goal in evolutionary biology, with important implications for conservation, medicine, and agriculture. Adaptation often leads to the fixation of beneficial alleles, and therefore it erodes local diversity so that understanding the coexistence of multiple adaptive phenotypes requires deciphering the ecological mechanisms that determine their respective benefits. Here, we show how antagonistic frequency-dependent selection (FDS), generated by natural and sexual selection acting on the same trait, maintains mimicry polymorphism in the toxic butterfly Heliconius numata Positive FDS imposed by predators on mimetic signals favors the fixation of the most abundant and best-protected wing-pattern morph, thereby limiting polymorphism. However, by using mate-choice experiments, we reveal disassortative mate preferences of the different wing-pattern morphs. The resulting negative FDS on wing-pattern alleles is consistent with the excess of heterozygote genotypes at the supergene locus controlling wing-pattern variation in natural populations of H. numata The combined effect of positive and negative FDS on visual signals is sufficient to maintain a diversity of morphs displaying accurate mimicry with other local prey, although some of the forms only provide moderate protection against predators. Our findings help understand how alternative adaptive phenotypes can be maintained within populations and emphasize the need to investigate interactions between selective pressures in other cases of puzzling adaptive polymorphism.

Evolutionary Novelty in a Butterfly Wing Pattern through Enhancer Shuffling.

An important goal in evolutionary biology is to understand the genetic changes underlying novel morphological structures. We investigated the origins of a complex wing pattern found among Amazonian Heliconius butterflies. Genome sequence data from 142 individuals across 17 species identified narrow regions associated with two distinct red colour pattern elements, dennis and ray. We hypothesise that these modules in non-coding sequence represent distinct cis-regulatory loci that control expression of the transcription factor optix, which in turn controls red pattern variation across Heliconius. Phylogenetic analysis of the two elements demonstrated that they have distinct evolutionary histories and that novel adaptive morphological variation was created by shuffling these cis-regulatory modules through recombination between divergent lineages. In addition, recombination of modules into different combinations within species further contributes to diversity. Analysis of the timing of diversification in these two regions supports the hypothesis of introgression moving regulatory modules between species, rather than shared ancestral variation. The dennis phenotype introgressed into Heliconius melpomene at about the same time that ray originated in this group, while ray introgressed back into H. elevatus much more recently. We show that shuffling of existing enhancer elements both within and between species provides a mechanism for rapid diversification and generation of novel morphological combinations during adaptive radiation.

Warning signals are under positive frequency-dependent selection in nature.

Positive frequency-dependent selection (FDS) is a selection regime where the fitness of a phenotype increases with its frequency, and it is thought to underlie important adaptive strategies resting on signaling and communication. However, whether and how positive FDS truly operates in nature remains unknown, which hampers our understanding of signal diversity. Here, we test for positive FDS operating on the warning color patterns of chemically defended butterflies forming multiple coexisting mimicry assemblages in the Amazon. Using malleable prey models placed in localities showing differences in the relative frequencies of warningly colored prey, we demonstrate that the efficiency of a warning signal increases steadily with its local frequency in the natural community, up to a threshold where protection stabilizes. The shape of this relationship is consistent with the direct effect of the local abundance of each warning signal on the corresponding avoidance knowledge of the local predator community. This relationship, which differs from purifying selection acting on each mimetic pattern, indicates that predator knowledge, integrated over the entire community, is saturated only for the most common warning signals. In contrast, among the well-established warning signals present in local prey assemblages, most are incompletely known to local predators and enjoy incomplete protection. This incomplete predator knowledge should generate strong benefits to life history traits that enhance warning efficiency by increasing the effective frequency of prey visible to predators. Strategies such as gregariousness or niche convergence between comimics may therefore readily evolve through their effects on predator knowledge and warning efficiency.

Mimicry among Unequally Defended Prey Should Be Mutualistic When Predators Sample Optimally.

Understanding the conditions under which moderately defended prey evolve to resemble better-defended prey and whether this mimicry is parasitic (quasi-Batesian) or mutualistic (Müllerian) is central to our understanding of warning signals. Models of predator learning generally predict quasi-Batesian relationships. However, predators' attack decisions are based not only on learning alone but also on the potential future rewards. We identify the optimal sampling strategy of predators capable of classifying prey into different profitability categories and contrast the implications of these rules for mimicry evolution with a classical Pavlovian model based on conditioning. In both cases, the presence of moderately unprofitable mimics causes an increase in overall consumption. However, in the case of the optimal sampling strategy, this increase in consumption is typically outweighed by the increase in overall density of prey sharing the model appearance (a dilution effect), causing a decrease in mortality. It suggests that if predators forage efficiently to maximize their long-term payoff, genuine quasi-Batesian mimicry should be rare, which may explain the scarcity of evidence for it in nature. Nevertheless, we show that when moderately defended mimics are profitable to attack by hungry predators, then they can be parasitic on their models, just as classical Batesian mimics are.

Genetic architecture and balancing selection: the life and death of differentiated variants.

Balancing selection describes any form of natural selection, which results in the persistence of multiple variants of a trait at intermediate frequencies within populations. By offering up a snapshot of multiple co-occurring functional variants and their interactions, systems under balancing selection can reveal the evolutionary mechanisms favouring the emergence and persistence of adaptive variation in natural populations. We here focus on the mechanisms by which several functional variants for a given trait can arise, a process typically requiring multiple epistatic mutations. We highlight how balancing selection can favour specific features in the genetic architecture and review the evolutionary and molecular mechanisms shaping this architecture. First, balancing selection affects the number of loci underlying differentiated traits and their respective effects. Control by one or few loci favours the persistence of differentiated functional variants by limiting intergenic recombination, or its impact, and may sometimes lead to the evolution of supergenes. Chromosomal rearrangements, particularly inversions, preventing adaptive combinations from being dissociated are increasingly being noted as features of such systems. Similarly, due to the frequency of heterozygotes maintained by balancing selection, dominance may be a key property of adaptive variants. High heterozygosity and limited recombination also influence associated genetic load, as linked recessive deleterious mutations may be sheltered. The capture of deleterious elements in a locus under balancing selection may reinforce polymorphism by further promoting heterozygotes. Finally, according to recent genomewide scans, balanced polymorphism might be more pervasive than generally thought. We stress the need for both functional and ecological studies to characterize the evolutionary mechanisms operating in these systems.

Chromosomal rearrangements maintain a polymorphic supergene controlling butterfly mimicry.

Supergenes are tight clusters of loci that facilitate the co-segregation of adaptive variation, providing integrated control of complex adaptive phenotypes. Polymorphic supergenes, in which specific combinations of traits are maintained within a single population, were first described for 'pin' and 'thrum' floral types in Primula and Fagopyrum, but classic examples are also found in insect mimicry and snail morphology. Understanding the evolutionary mechanisms that generate these co-adapted gene sets, as well as the mode of limiting the production of unfit recombinant forms, remains a substantial challenge. Here we show that individual wing-pattern morphs in the polymorphic mimetic butterfly Heliconius numata are associated with different genomic rearrangements at the supergene locus P. These rearrangements tighten the genetic linkage between at least two colour-pattern loci that are known to recombine in closely related species, with complete suppression of recombination being observed in experimental crosses across a 400-kilobase interval containing at least 18 genes. In natural populations, notable patterns of linkage disequilibrium (LD) are observed across the entire P region. The resulting divergent haplotype clades and inversion breakpoints are found in complete association with wing-pattern morphs. Our results indicate that allelic combinations at known wing-patterning loci have become locked together in a polymorphic rearrangement at the P locus, forming a supergene that acts as a simple switch between complex adaptive phenotypes found in sympatry. These findings highlight how genomic rearrangements can have a central role in the coexistence of adaptive phenotypes involving several genes acting in concert, by locally limiting recombination and gene flow.

Into the Andes: multiple independent colonizations drive montane diversity in the Neotropical clearwing butterflies Godyridina.

Understanding why species richness peaks along the Andes is a fundamental question in the study of Neotropical biodiversity. Several biogeographic and diversification scenarios have been proposed in the literature, but there is confusion about the processes underlying each scenario, and assessing their relative contribution is not straightforward. Here, we propose to refine these scenarios into a framework which evaluates four evolutionary mechanisms: higher speciation rate in the Andes, lower extinction rates in the Andes, older colonization times and higher colonization rates of the Andes from adjacent areas. We apply this framework to a species-rich subtribe of Neotropical butterflies whose diversity peaks in the Andes, the Godyridina (Nymphalidae: Ithomiini). We generated a time-calibrated phylogeny of the Godyridina and fitted time-dependent diversification models. Using trait-dependent diversification models and ancestral state reconstruction methods we then compared different biogeographic scenarios. We found strong evidence that the rates of colonization into the Andes were higher than the other way round. Those colonizations and the subsequent local diversification at equal rates in the Andes and in non-Andean regions mechanically increased the species richness of Andean regions compared to that of non-Andean regions ('species-attractor' hypothesis). We also found support for increasing speciation rates associated with Andean lineages. Our work highlights the importance of the Andean slopes in repeatedly attracting non-Andean lineages, most likely as a result of the diversity of habitats and/or host plants. Applying this analytical framework to other clades will bring important insights into the evolutionary mechanisms underlying the most species-rich biodiversity hotspot on the planet.

Refining mimicry: phenotypic variation tracks the local optimum.

Müllerian mimicry between chemically defended preys is a textbook example of natural selection favouring phenotypic convergence onto a shared warning signal. Studies of mimicry have concentrated on deciphering the ecological and genetic underpinnings of dramatic switches in mimicry association, producing a well-known mosaic distribution of mimicry patterns across geography. However, little is known about the accuracy of resemblance between natural comimics when the local phenotypic optimum varies. In this study, using analyses of wing shape, pattern and hue, we quantify multimodal phenotypic similarity between butterfly comimics sharing the so-called postman pattern in different localities with varying species composition. We show that subtle but consistent variation between populations of the localized species, Heliconius timareta thelxinoe, enhance resemblance to the abundant comimics which drive the mimicry in each locality. Those results suggest that rarer comimics track the changes in the phenotypic optimum caused by gradual changes in the composition of the mimicry community, providing insights into the process by which intraspecific diversity of mimetic pattern may arise. Furthermore, our results suggest a multimodal evolution of similarity, with coordinated convergence in different features of the phenotype such as wing outline, pattern and hue. Finally, multilocus genotyping allows estimating local hybridization rates between H. timareta and comimic H. melpomene in different populations, raising the hypothesis that mimicry refinement between closely related comimics may be enhanced by adaptive introgression at loci modifying the accuracy of resemblance.

Selection of Valid Reference Genes for Reverse Transcription Quantitative PCR Analysis in Heliconius numata (Lepidoptera: Nymphalidae).

Identifying the genetic basis of adaptive variation is challenging in non-model organisms and quantitative real time PCR. is a useful tool for validating predictions regarding the expression of candidate genes. However, comparing expression levels in different conditions requires rigorous experimental design and statistical analyses. Here, we focused on the neotropical passion-vine butterflies Heliconius, non-model species studied in evolutionary biology for their adaptive variation in wing color patterns involved in mimicry and in the signaling of their toxicity to predators. We aimed at selecting stable reference genes to be used for normalization of gene expression data in RT-qPCR analyses from developing wing discs according to the minimal guidelines described in Minimum Information for publication of Quantitative Real-Time PCR Experiments (MIQE). To design internal RT-qPCR controls, we studied the stability of expression of nine candidate reference genes (actin, annexin, eF1α, FK506BP, PolyABP, PolyUBQ, RpL3, RPS3A, and tubulin) at two developmental stages (prepupal and pupal) using three widely used programs (GeNorm, NormFinder and BestKeeper). Results showed that, despite differences in statistical methods, genes RpL3, eF1α, polyABP, and annexin were stably expressed in wing discs in late larval and pupal stages of Heliconius numata This combination of genes may be used as a reference for a reliable study of differential expression in wings for instance for genes involved in important phenotypic variation, such as wing color pattern variation. Through this example, we provide general useful technical recommendations as well as relevant statistical strategies for evolutionary biologists aiming to identify candidate-genes involved adaptive variation in non-model organisms.