Perfect mimicry between Heliconius butterflies is constrained by genetics and development

Müllerian mimicry strongly exemplifies the power of natural selection. However, the exact measure of such adaptive phenotypic convergence and the possible causes of its imperfection often remain unidentified. Here, we first quantify wing colour pattern differences in the forewing region of 14 co-mimetic colour pattern morphs of the butterfly species Heliconius erato and Heliconius melpomene and measure the extent to which mimicking colour pattern morphs are not perfectly identical. Next, using gene-editing CRISPR/Cas9 KO experiments of the gene WntA, which has been mapped to colour pattern diversity in these butterflies, we explore the exact areas of the wings in which WntA affects colour pattern formation differently in H. erato and H. melpomene. We find that, while the relative size of the forewing pattern is generally nearly identical between co-mimics, the CRISPR/Cas9 KO results highlight divergent boundaries in the wing that prevent the co-mimics from achieving perfect mimicry. We suggest that this mismatch may be explained by divergence in the gene regulatory network that defines wing colour patterning in both species, thus constraining morphological evolution even between closely related species.


Introduction
Adaptation is the product of natural selection on the genetic and phenotypic diversity within a population [1]. In this regard, key developmental steps that limit or bias trait variation can pose so-called constraints on the directionality of evolution [2][3][4]. Consequently, when populations evolve independently, lineages can accumulate changes that lead evolution along an irreversible trajectory [5]. Understanding the relative contribution (or constraint) of genetics and development to adaptation would therefore allow us to better comprehend the directionality and predictability of evolution [6].
To date, the extent to which evolutionary changes are consequential for future adaptation has been most elegantly studied using artificial selection experiments [1]. In these experiments, multiple generations can be relatively easily traced while being exposed to contrasting selection pressures. The consequences of their adaptations can then be investigated when these selection pressures are reversed. In malaria, for example, a single mutation of large effect in the K13 protein can confer drug resistance [7] but simultaneously also favours the evolution of additional epistatic mutations [8]. As a consequence of these changes, the evolved phenotype cannot be simply reversed to its ancestral state when withdrawing anti-malarial medication and thus has important consequences for resistance management strategies [8]. In butterflies, artificial selection on eyespots has been used to test for the existence of developmental constraints in wing colour patterns. Because these wing pattern elements belong to the same homologous series and share developmental pathways, covariances between them are expected in their size and shape. Nevertheless, in the butterfly Bicyclus anynana, artificial selection has suggested that there is great potential for independent change in size and shape of different eyespots [9]. This observation has been used to argue that natural selection plays a dominant role over developmental constraints in the observed eyespot diversity among Bicyclus species [9]. Similarly, in Drosophila, artificial selection experiments have shown that changes to wing shape can be rapidly induced from standing genetic variation present in the populations [10]. However, these induced phenotypes are lost when selection is suspended, presumably due to pleiotropic links to other traits that result in negative fitness consequences [11].
Not all study systems are amenable to perform artificial selection experiments. Alternatives to such controlled experiments are cases of convergent evolution that provide natural opportunities to investigate the selective, genetic, and developmental routes to adaptation [1]. For example, mimicry and the resulting phenotypic convergent evolution between distinct butterfly species provides a comparative framework to investigate the genes underlying the evolution and diversity of a wing colour pattern [12][13][14]. Recent studies have shown, for example, that the genes WntA, cortex, and optix are repeatedly used to control variation in wing colour patterns across Nymphalid butterflies [14][15][16][17]. Nevertheless, even in cases of Müllerian mimicry between species within the Heliconius genus, in which both partners have used the same genes to converge on an aposematic warning signal, some degree of imperfection in resemblance may exist. However, the precise extent to which Heliconius mimetic butterflies need to perfectly resemble the same phenotype to maximize the fitness value is not well understood. What may underlie these differences in resemblance are (i) conflicting or relaxed selection pressures and/or (ii) genetic and developmental constraints. Conflicting selective pressures can include variation in the mimicry community [18] and conflict between the outcomes of natural and sexual selection [19]. Relaxed selection pressures may result from coarse discrimination by predators [20,21]. On the other hand, genetic and developmental constraints can result from divergence in the genetic background or in the assembled gene regulatory network that affects the functioning and phenotypic effect of these genes.
Heliconius erato and Heliconius melpomene provide a textbook case of Müllerian mimicry. Although there is no evidence for gene flow between H. erato and H. melpomene, which split around 12-14 Mya [22], their resemblance in wing colour patterns is remarkable. This phenotypic convergence has evolved through strong selection pressures that benefit a common warning pattern that birds have learned to associate with unpalatability [23,24]. The discriminatory visual properties of birds appear to be quite precise, resulting in strong selection pressures for fine scale adjustments of the shape and size of colour patterns among local mimetic butterfly communities (for an overview of selection coefficients see [25]) [26,27].
Recently, a series of functional experiments have tested the role of the WntA gene in different Heliconius butterfly species and populations [12]. WntA is a member of the Wnt family of signalling ligands and a key molecular tool for butterfly wing colour pattern development. While the gene coding sequence is highly conserved across Lepidoptera, its cis-regulatory diversity underlies wing pattern shape variation between and within butterflies species [28,29]. Recent CRISPR/Cas9 KO experiments have shown its role in defining the ultimate colour fate of individual wing scale cells and highlighted incredible variability in the position and wing territory affected by WntA in H. erato and H. melpomene [12]. This result provides a visual representation of the effect of a divergent genetic background on a gene's function. In H. erato and H. melpomene, CRISPR/Cas9 WntA mutant phenotypes evidenced a more restricted area of modified black scales in the forewing of H. melpomene compared to H. erato [12,17]. This result demonstrates that although the two co-mimetic butterflies display identical forewing colourations, the genetic architecture and gene regulatory networks underlying their resemblance may be more complex than previously thought. These differences likely arise from divergence in the regulation of WntA and/or divergence in epistatic interactions with other genes, together defined as the gene regulatory network [12].
In the current study, we aimed to explore the contribution of genetic and developmental constraints to the phenotypic mismatch between co-mimetic pairs of Heliconius butterflies. We therefore first quantitatively measured and compared differences in the mid-forewing band (MFB) pattern of 14 co-mimicking populations of the butterfly species H. erato and H. melpomene from Central and South America. We then interpreted these differences in light of the recent Heliconius WntA CRISPR/Cas9 KO mutants [12]. We were specifically interested to test how species-specific divergence might impact the developmental function of WntA and limit adaptive convergence. More precisely, we argue that overlap of wild-type differences with differences in pattern boundaries as defined by WntA KOs in H. erato and H. melpomene would suggest the possible existence of genetic and developmental constraints for natural selection to achieve perfect mimicry. By combining the most recent methodological advances in functional experiments with quantitative measures, our study offers an alternative approach to artificial selection experiments to test the relative constraint of genetics and development to adaptation.

Materials (a) Colour pattern analysis
We obtained 8-14 images of each of 14 mimicking colour pattern morphs (electronic supplementary material, tables S1 and S2). Images were obtained through the authors' collections, collaborations, and collections made publicly available by Cuthill et al. [30] and Jiggins et al. [31]. Individual genders were determined based on sexual dimorphism in the androconial region [32]. To align images, we used a total of 11 landmarks at vein intersections on one forewing of each sample (electronic supplementary material, figure S1). Methods and results describing additional analysis to control for interspecific changes in wing shape and sex differences are described in Supplementary materials S1 and electronic supplementary material, table S3. While WntA is involved in black scale development in both the fore-and hindwing in Heliconius [12,17], differences in the distribution of WntA have mainly been found to correlate with the position of black royalsocietypublishing.org/journal/rspb Proc. R. Soc. B 287: 20201267 colour in the central part of the forewing among Heliconius colour pattern morphs, here called mid-forewing band or MFB [28]. With the interest of studying variation in the MFB, we extracted and focused on the area of the forewing in which black is absent and in which WntA is thus likely not expressed. MFB patterns were extracted and aligned using the R package patternize [33]. Depending on the MFB phenotype, we specified red, green, and blue (RGB) values for red, yellow, and/or white with a colour threshold (colOffset) chosen to fully extract the pattern. For H. e. notabilis, H. m. plesseni, and H. m. cythera, we extracted and combined both red and white to represent the MFB shape. Background noise or damaged regions in the wing that were co-extracted with the colour patterns were masked using the setMask function. Next, a thin plate spline (tps) transformation was obtained from transforming landmarks to a common reference sample. This common reference sample included the landmarks of an arbitrarily chosen sample and was used in all colour pattern analysis. The tps transformation was then used to align and compare the extracted MFB shape, size, and position.
Differences in the MFB patterns were first compared by subtracting the H. erato and H. melpomene pattern frequencies of each population, obtained with the sumRaster function in patternize (i.e. absolute MFB difference) and compared between co-mimics using a one-sample t-test. Next, the relative size of the pattern was calculated as the proportion of the total wing area in which the pattern is observed, using the patArea function (i.e. relative MFB difference) and compared between co-mimics using a two-sample t-test. Principal component analysis (PCA) was performed on the binary representation of the aligned colour pattern rasters of each sample [33]. The PCA visualizes the main variations in colour pattern boundaries among samples and groups and provides predictions of colour pattern changes along the PC axis, with positive values presenting a higher predicted presence of the MFB pattern and negative values presenting the absence of the pattern. Parts of the colour patterns that are present in all samples have a predicted value of zero, as these pixels do not contribute variance in the PCA. We tested the effect of population, sex, and species on shape variables using multivariate analysis of variance (MANOVA) and linear discriminant analysis (LDA) using only the values of samples along significant PC axes (see Supplementary Materials S1 for details).

(b) WntA CRISPR KO analysis
Five mutant butterflies (10 wings) for each of the Panamanian geographical colour pattern morphs Heliconius erato demophoon and Heliconius melpomene rosina for which a frame shift mutation was generated at the gene WntA using CRISPR/Cas9 were obtained from Concha et al. [12]. All these mutants showed symmetric changes in wing patterns on both the left and right forewings and were thus likely full KO mutants [12]. Both left and right forewings were landmarked and the mutant pattern was extracted and aligned using the R package patternize [33]. Red was extracted from the H. e. demophoon mutants. As the H. m. rosina mutants often showed a yellow spot appearing in the proximal part of the MFB, both red and yellow were extracted for H. m. rosina. The mutant patterns were superimposed on the wild-type wing pattern comparisons by aligning to the common reference sample and using the contour function of the R package raster [34].

Results
(a) Divergence and convergence in mid-forewing band pattern Geographical colour pattern morphs of H. erato and H. melpomene cover a wide spectrum of MFB patterns, with unique or partially overlapping pattern elements among them (figure 1a,b). In the PCA of the MFB, the first main axis of variance was dominated by the absence or presence of a broad red MFB, also typically called a 'Postman' phenotype (PC1, 32% of variation; figure 1c).  figure 1d). This difference between species was highly significant along PC2 (10% of variation; F 1,127 = 131.1, p < 0.001) and suggests a general trend of expansion of both the proximal and distal area of the MFB in H. melpomene Postman phenotypes compared to H. erato. Significant differences between male and female MFB patterns were observed in both H. erato (F 1,139 = 3.11, p = 0.002) and H. melpomene (F 1,140 = 3.17, p = 0.001). However, sex differences had a low probability of posterior classification (55.8% and 58.3% for male and female H. erato and 78.0% and 62.5% for male and female H. melpomene, which is close to classification between random groups) and were only significant along PC axes that explain small amounts of variation among samples (electronic supplementary material, table S4).
Using pairwise comparisons we further explored interspecific differences between H. erato and H. melpomene co-mimics

(b) Mismatch between co-mimics partly coincides with developmental WntA boundaries
An important aspect of our study was to determine the possible link between wing pattern variation observed between wild-type butterflies and the WntA CRISPR/Cas9 KO phenotypes. We focused on the Panamanian co-mimics H. e. demophoon and H. m. rosina for which the largest WntA KO dataset is available and compared the WntA boundaries defined by these mutants with wild-type variation in the Postman phenotypes. As recently described by Concha et al.

Discussion
In Heliconius butterflies, convergence between co-mimicking populations has been broadly defined as nearly identical. In accordance with these phenotypic similarities, genetic work has suggested that convergence is governed by a small and shared set of genes [13,15,17,28,36]. For example, WntA has been linked to wing colour pattern variation in the forewing of both H. erato and H. melpomene butterflies [28]. Recently, the availability and implementation of CRISPR/Cas9 functional approaches in butterfly wing colour pattern development has allowed to further investigate the role of wing patterning genes [12,17,37]  royalsocietypublishing.org/journal/rspb Proc. R. Soc. B 287: 20201267 developmental background. However, the extent to which divergence in WntA functioning may have affected the evolution of natural populations of mimetic phenotypes in Heliconius has not been investigated. By doing so, our study provides new insights into the dynamics of selection and how adaptation may be constrained to find alternative routes.

(a) Patterns of advergence in mid-forewing band
A question that arises from the observed differences in MFB patterns between co-mimics is whether it is necessary for H. erato and H. melpomene to perfectly mimic each other and whether selective forces could work against convergence [19,38]. In our comparisons we observed that even though the position of forewing pattern elements may not be perfectly identical between co-mimics, the relative amount of black versus red or yellow is generally more similar than the absolute difference. This improved match of the size of the MFB seems to result from compensatory changes in the MFB pattern and are in line with phenotypic evolution being driven by predation pressure.  [27]. They thus indicate fine scale pattern adaptation in non-homologous regions of the wing to obtain a better match in the shape and size of the pattern even though they have a shifted position in the wings. Differences in MFB patterns between sexes were the same for both species and suggest that sexual conflicts for species recognition are likely not driving species differences in MFB patterns. The compensatory evolution to obtain a more similar area of the MFB despite its mismatch in position may suggest imperfect discrimination in the visual range of their predators, or the relative importance of overall features of colour contrast distribution rather than the exact position of pattern elements [20,21]. A remarkable example of this are the co-mimetic Ecuadorian butterflies H. e. notabilis and H. m. plesseni which both have red and white in their proximal forewing element but have the relative positions of white versus red colour inversed. Notably, these wing colour patterns are the results of complex epistatic interactions between WntA and other genes such as the transcription factor optix, which controls white and red scale development [13,36].
We also found detectable within-species differences between populations that are generally considered identical. Among Postman populations, the greatest MFB differences were observed between several H. erato populations (e.g. between H. e. hydara from Panama and French Guiana). Notably, these H. erato Postman populations also had the most apparent differences with its H. melpomene co-mimic (i.e. H. e. hydara versus H. m. melpomene from French Guiana). As H. erato is often suggested to be the more abundant co-mimic and, thus, the model which H. melpomene mimics [21,22], the presence of a less perfect mimetic signal in these cases may reflect a 'lag' in the evolution of better mimicry of the H. melpomene populations. This inference potentially fits a signal of recent adaptive evolution (i.e. selective sweep signal) across the regulatory regions in the first intronic region of WntA of H. e. hydara populations from French Guiana, but not H. e. hydara populations from Panama, or H. m. melpomene populations from French Guiana [25]. What drives the divergence in MFB pattern between the H. erato Postman populations may include local changes in the composition of the mimicry and/or predator community [18,24,26].

(b) Indications of genetic and developmental constraints
Some of the pattern differences we observed between comimetic populations were shared among several geographically distinct populations and phenotypes. Such shared differences thus likely result from a shared genetic and developmental background of specific MFB areas among these geographical populations. This observation highlights the existence of a WntA KO boundary and suggests that the observed imperfections in mimicry might be to some extent imposed by divergence in the gene regulatory network involved in the development of the MFB. As discussed earlier, it is possible that conflicting selection pressures can explain imperfections in mimicry. However, our quantitative wing colour pattern analyses on wild-type individuals and WntA KO phenotypes suggests the likely role of divergent genetic backgrounds in constraining the observed imperfect mimicry. While not obvious in the pairwise comparisons of the Postman co-mimics, the PCA analysis identified that the distal anterior margin of the MFB is generally expanded in H. melpomene compared to H. erato. Finally, the difference observed in the distal anterior margin of the MFB between co-mimicking Postman populations matched the distal boundary of WntA affected wing area as identified in the H. e. demophoon WntA KOs (figure 3c) and may potentially further highlight areas of the wing that are constrained by divergence in the gene regulatory network between H. erato and H. melpomene.

(c) Candidates of divergence in the gene regulatory network
Divergence in the gene regulatory network that is involved in the development of the MFB may include a multitude of changes in upstream factors that regulate spatial and/or temporal expression of WntA [39], the cis-regulatory elements of WntA itself [29], as well as additional genes that define the fate of wing scale cells. A few of these diverged elements are likely loci or genes that have previously been implicated in  [40] and these loci have been demonstrated to include cis-regulatory elements of WntA [29]. In H. melpomene, on the other hand, regulatory variation at the so-called N locus, which includes the gene cortex and a few additional candidate genes, seems to underlie MFB shape variation together with WntA [40]. Potentially, this latter locus provides a candidate that explains the absence of a WntA effect on black wing scale development in the proximal area of the wing in H. melpomene. Finally, apart from the major effect loci involved in Heliconius colour patterns, quantitative trait loci (QTL) studies of pattern variation in H. erato and H. melpomene have demonstrated the existence of minor effect loci associated with quantitative changes in wing colour pattern [42][43][44]. This larger set of genetic variants controlling quantitative variation is additional to the regulatory complexity that modulates the expression of the major colour pattern genes [12,29,37].

Conclusion
Studying the evolutionary dynamics of adaptation is complicated due to the interplay of selection, genetics, and development. Adaptive radiations and mimicry systems have long been powerful study systems to unravel the genetic basis of adaptive traits. Here, we used Heliconius butterflies to investigate the interplay of selection and genetics and provide a tentative explanation on the causes that limit perfect wing mimicry after several million years of strong natural selection towards a mutual anti-predatory signal. We propose that phenotypic patterns of adaptation and convergence between Heliconius co-mimics is biased by divergence in the gene regulatory network underlying the mimicry trait. These networks can be highly polygenic and interconnected and evolution may not have been able to perfectly rewire them in the same identical way after the accumulation of changes across their genomes.
In the case of the observed differences between the WntA mutant phenotypes of H. erato and H. melpomene, these constraints may exist due to independently evolved genetic elements that interact with the WntA gene or protein. These elements may include trans factors, cis-regulatory elements, or additional genes with a complementary role to WntA and may not have the genetic architecture to easily be detached from potential developmental interactions with other genes [45,46]. While we quantified exact wing colour pattern differences between H. erato and H. melpomene co-mimics, we also observed that the warning signal represented by the relative amount of colours is more consistent. Hence, changes to the MFB can be observed that compensate the developmental bias in other areas of the wing, indicating that selection has used an alternative route to convergence. Overall, our results provide a quantitative measure of the imperfection of mimicry and propose a novel view into the interplay between selective conflicts compared to genetic and developmental constraints in the production of similar phenotypes.