Convergent morphology and divergent phenology promote the coexistence of Morpho butterfly species

The coexistence of closely-related species in sympatry is puzzling because ecological niche proximity imposes strong competition and reproductive interference. A striking example is the widespread wing pattern convergence of several blue-banded Morpho butterfly species with overlapping ranges of distribution. Here we perform a series of field experiments using flying Morpho dummies placed in a natural habitat. We show that similarity in wing colour pattern indeed leads to interspecific territoriality and courtship among sympatric species. In spite of such behavioural interference, demographic inference from genomic data shows that sympatric closely-related Morpho species are genetically isolated. Mark-recapture experiments in the two most closely-related species unravel a strong temporal segregation in patrolling activity of males. Such divergence in phenology reduces the costs of reproductive interference while simultaneously preserving the benefits of convergence in non-reproductive traits in response to common ecological pressures. Henceforth, the evolution of multiple traits may favour species diversification in sympatry by partitioning niche in different dimensions.

N atural communities are composed of multiple species involved in diverse ecological interactions either facilitating or impairing their coexistence. These interactions strongly depend on the level of phylogenetic divergence between the involved species 1,2 . In particular, the coexistence of closelyrelated species is often impaired by their ecological and morphological similarities [3][4][5] , generating either elevated interspecific competition and/or reproductive interference [6][7][8] . Close relatedness and shared phenotypes are indeed identified as important factors enhancing behavioural interference in sympatric animal species 9 . Nevertheless, inherited traits shared by closely-related species may also provide individual fitness benefits, for instance when facing common predators either by sharing a common warning signal 10,11 , or by responding to alarm cues emitted by heterospecifics 12,13 . Such interspecific positive densitydependence is predicted to favour species coexistence 14 . Here, we investigate how shared and divergent traits can facilitate the coexistence of closely-related species.
Mimetic butterflies are a striking case of sympatric species with phenotypic convergences strongly promoted by natural selection. The convergence observed for warning wing patterns in defended species is frequently followed by convergence in other traits like flight height 15,16 or host-plant 17 since sharing a microhabitat increases the similarity of encountered predatory communities, thereby enhancing protection 18 . Nevertheless, the benefits conferred by overlapping visual signals and ecological niches in mimetic species may, in turn, incur fitness costs through increased heterospecific rivalry, heterospecific female harassment and the expression of Dobzhansky-Müller incompatibilities in hybrids [19][20][21] . These costs might ultimately either limit the coexistence of closely-related mimetic species 22 or promote the evolution of alternative cues involved in species recognition 23 . The conflicting ecological interactions between mimetic species, therefore, question their persistence in sympatry when they are closely related and point to the evolution of alternative divergent traits.
Here we investigate how trait evolution might favour the coexistence of closely-related species by focusing on three species of the butterfly genus Morpho that are sympatric over most of the Amazonia 24,25 and occupy the same understory niche in the tropical forest. Multiple local convergences in wing colour patterns are observed among these species, throughout their geographical range 25 . Although not chemically defended, these butterflies are very difficult to capture because of their fast, erratic flight. The contrast between their dorsal bright iridescent blue and ventral cryptic brownish wing surfaces induces a flash pattern during flapping flight, that was suggested to confuse predators, further increasing their difficulty of capture 26 . Predators may then learn to avoid such elusive prey harbouring the conspicuous blue patterns. The local convergence in the three closely-related Morpho species probably stems from frequency-dependent selection generated by predator behaviour, in a similar way as in Müllerian mimics 25 . The iridescent blue colouration of Morpho butterflies shared by sympatric species may thus reduce individual predation by advertising escape ability (the 'escape mimicry hypothesis', see ref. 27 ). Males and females typically display the same colour pattern shared between sympatric species and may thus benefit from increased protection against the same predator community. While females are rarely observed in the field, males from these different closely-related species typically patrol within the same habitats. Male-male interactions are then frequently observed in these butterflies, one of the males eventually being chased away. Male territoriality may therefore occur both within and among species, potentially leading to important reproductive interferences 28 . The local convergence in colour pattern might further enhance heterospecific rivalry and impair mate recognition in these sympatric species. The persistence of closelyrelated species sharing similar colouration would thus be enabled if intraspecific benefits of maintaining a convergent colour pattern outweigh the costs associated with reproductive interference.
To test whether convergent wing patterns generate reproductive interference between species, we investigated the behaviour of Morpho towards flying dummy butterflies harbouring various wing patterns, in a series of experiments performed in the wild. To infer whether these heterospecific reproductive behaviours underpin genetic exchange, we then estimated the gene flow among them, their level of genomic divergence and investigated their demographic history. To explore whether interference could be mitigated by temporal partitioning, we finally monitored male flight activity using capture-mark-recapture experiments.
Our study reveals strong genetic isolation between species despite reproductive interference, probably facilitated by a marked difference in flying hours between Morpho species. This work highlights how joint convergent and divergent evolution of different traits may promote species diversification in sympatry.

Results and discussion
We focused on a single locality from Amazonian Peru where three Morpho species with strikingly similar colour patterns (Morpho achilles, Morpho helenor and Morpho deidamia), live in sympatry. In this field site, males from these three closely-related species display a typical patrolling behaviour along the river bed, allowing us to investigate species interactions in natura.
Heterospecific interactions lead to reproductive interference among sympatric species. To test whether the strong convergence in colour pattern among these sympatric Morpho species leads to heterospecific rivalry and courtship, we investigated the response of patrolling males to butterfly dummies placed in the field. We built realistic moving butterfly dummies using the actual wings of captured Morpho, set up on a solar-powered fluttering device ( Fig. 1a and Supplementary Video 1; see also Supplementary Figs. 1, 2). To ensure that only visual cues were triggering the interactions, the wings were washed in hexane prior mounting, ruling out any effect of pheromones 29 . We built ten different dummies with the wings of specimens from different species and sexes. Eight dummies were built from butterflies caught in the same Peruvian site, using both sexes of four sympatric species: the three mimetic species M. helenor, its sister species, M. achilles, and the more distantly-related M. deidamia, all exhibiting an iridescent blue band bordered by proximal and distal black areas, and the phenotypically distinct Morpho menelaus, exhibiting fully blue iridescent wings (referred to as local dummies: n = 8) (Fig. 1b, c and Supplementary Fig. 1). The last two dummies were built from M. helenor and M. achilles males captured in French Guiana (referred to as exotic dummies: n = 2), exhibiting a narrower blue band relative to local-Peruvian-individuals ( Fig. 1b and Supplementary Fig. 1). We tested all dummies (n = 10) following a randomised design: a different dummy was placed each morning at the same site and left fluttering on the river bank for 5 h. This was replicated four times per dummy. The dummy was continuously filmed using a camera (Gopro Hero5 Black set at 120 images per second) and monitored by a human observer who recorded the timing of responses displayed by wild butterflies.
Over a 2-months period, we recorded the patrolling behaviours of all males passing through the river on sunny mornings, resulting in 2700 responses to dummies. We specifically focused on the behaviour of butterflies from the two mimetic sister species M. helenor and M. achilles, which represented the large majority of the passing males (35% and 47% respectively, Supplementary  Fig. 3). Visiting species were identified by sight, except for M. helenor and M. achilles that cannot be reliably discriminated in flight. We used an indirect method to distinguish these two strikingly similar species (see 'Methods' section and Supplementary Fig. 13). During these sessions, we defined two behaviours: (1) approaches, when a marked change in the trajectory of the passing butterfly towards the dummy was observed and (2) interactions, when the butterfly entered a 1 m 3 zone around the dummy (Fig. 1a). We controlled for the effect of cloud cover on these behaviours (see 'Methods' section and Supplementary  Fig. 4). This procedure allowed us to test whether patrolling Morpho males (1) are more strongly attracted by the colour pattern of their conspecifics as compared to that of other species, (2) discriminate between sexes and (3) are more attracted by local colour patterns than by exotic ones, thereby assessing the strength of behavioural interference between sympatric species.
About half of the patrolling individuals deviated from their flight path to approach the setup (see Supplementary Figs Table 7): the long-range blue signal emitted by the different wing patterns displayed in the different species thus appeared similarly attractive to patrolling M. helenor and M. achilles males, suggesting they may approach anything roughly recognised as a potential mate or rival.
In both M. helenor and M. achilles, about 40% of approaches resulted in interactions, the visitor typically flying in circles around the dummy (Supplementary Video S1). Surprisingly, interactions with conspecific dummies were not significantly higher than with sympatric heterospecifics. Overall, no strong discrimination between sympatric conspecifics and sympatric heterospecifics, nor between female and male dummies was found, suggesting that interspecific interactions are frequent among wild Morphos in sympatry (Fig. 1b, c). Some differences were nevertheless detected. Dummies with a larger wing area and a greater proportion of iridescent blue colouration were more likely to trigger interactions ( Supplementary Fig. 7), suggesting that these characteristics are used by patrolling males to discriminate encountered individuals. M. helenor males were found to approach dummies (all sexes and species pooled) significantly more than M. achilles males (mean % of approach in M. helenor = 48.3 ± 15.7; M. achilles = 39.7 ± 14.8; P = 0.01). They interacted as frequently with M. achilles females as with their conspecific females (Fig. 1b). In contrast, interactions occurred in markedly lower proportions with the exotic dummies, that were largely ignored (Fig. 1b). Morpho males, therefore, do distinguish different colour patterns (i.e. large blue band in local individuals, vs. narrower blue band in exotic ones), but yet they largely engage in interactions with congeners bearing locally known signals, including signals sharply dissimilar from their own displayed by more distantly-related species (e.g. fully blue wings in M. menelaus). For males, interspecific interactions with both males and females could be promoted because mating opportunities are limited by strong male-male competition. Defending territory against any males harassing a potential mate, or courting any females would then be a favoured behaviour. The costs of missed mating opportunities are probably elevated for males when species recognition is impaired by the phenotypic resemblance between species 10 .
Overall, the strong visual attraction and the low species discrimination suggest that both heterospecific interactionscontests among males and mating or mating attempts with females-can frequently happen. This indicates that reproductive interferences occur among these sympatric Morpho species.
Patrolling M. achilles males had limited interactions with the co-mimetic M. helenor female, but interacted indiscriminately with the dummies of the other sympatric species (i.e. M. deidamia   c Frequency of interaction with the female dummies (yellow colour). Raw data << nb of interactions/nb of approaches >> are indicated on each bar. Proportions were compared using Fisher Exact probability tests. Only significant differences are shown (P < 0.05). For proportions differing from all the others, the highest P-value is reported. Source data are provided as a Source Data file.
and M. menelaus) (Fig. 1c). Such a specific, acute visual discrimination towards its -albeit phenotypically closest-sister species possibly evolved in M. achilles as a result of reinforcement selection against genetic incompatibilities expressed in hybrids 30 .
Whether reinforcement selection does occur in those closelyrelated species would however require further study. Reinforcement process may also occur through divergence in olfactory cues enabling discrimination among species 20,31 , although this remains to be investigated in Morpho. This specific effect detected in female interactions suggests that patrolling males can discriminate sexes based on short-distance visual cues. To test whether the behaviour of M. achilles males indeed differs when interacting with conspecific male and female dummies, we equipped our set up with a stereoscopic high-speed videography system, enabling us to quantify the threedimensional flight kinematics of visiting butterflies in natural conditions on a sub-set of sessions (Figs. 1a and 2a). Behavioural differences were detected between the interactions toward males and females (n = 14 flights analysed for each sex): wild males circled closer to the female compared to the male dummy (significant interaction between distance from dummy and the sex of the dummy on the proportion of time spent within distance intervals: F 11,312 = 1.9; P = 0.03, Fig. 2b, see also Supplementary  Fig. 8), ensuring a steady speed when close to the female, whereas they displayed more erratic accelerations around the male dummy (significant interaction between the distance from dummy and the sex of the dummy on acceleration within distance intervals: F 1,177 = 1.9; P < 0.001, Fig. 2c). These aerial interactions suggest that M. achilles males adjust their flight behaviour according to the sex of their conspecific, showing that, at least in M. achilles, Morpho males can identify females using visual cues alone.
The limited attraction of M. achilles towards the females of its sister species M. helenor and the contrasted attraction to local and exotic colour patterns suggest that sympatric interactions might have promoted discrimination capacities in the closest related species, thereby reinforcing pre-zygotic barriers. Overall, the behavioural interferences among these sympatric Morpho species as well as the putative reinforcement effect detected in M. achilles question their stable coexistence and reproductive isolation.
Limited genetic exchanges between species despite close relatedness and reproductive interference. To test the level of reproductive isolation despite such a strong behavioural interference, we then used genomic data obtained for the three mimetic species M. achilles, M. helenor and M. deidamia. Patterns of genetic polymorphism (Supplementary Figs. 9 and 10) and divergence ( Supplementary Fig. 11) throughout the genome were investigated using RAD-sequencing. Eight categories of contrasted scenarios of speciation and models of linked selection ( Supplementary Fig. 12) were then proposed to be statistically evaluated using a modified version of DILS 1.0.0 adapted to three gene pools 32 . Each category was simulated using four sub-models depending on whether the effective population size is assumed to be homogeneous or heterogeneous along the genome, as well as for introgression rates. Consistent with the phylogeny 33 , demographic inferences performed on our Peruvian populations revealed a more recent time of the split between the sister species M. achilles and M. helenor than between these two species and M. deidamia. Our hierarchical approach of eight categories of models according to temporal patterns of migration between all three species provided strong statistical support for current genetic isolation between M. deidamia and both M. helenor/M. achilles (F st ± std M. helenor -M. achilles = 0.30 ± 0.226; M. helenor -M. deidamia = 0.82 ± 0.11; M. achilles -M. deidamia = 0.85 ± 0.11) ( Fig. 3; posterior probability = 0.94). Current isolation was also strongly supported between M. helenor and M. achilles, with gene flow restricted to the early times of speciation (posterior probability = 0.80). Genomic or behavioural barriers are thus likely to strongly limit gene flow between these closely-related sympatric species (Fig. 3 and Supplementary Tables 8-10). Note that if our demographic inferences show that models with ancestral migration reproduce the molecular data better than models with ongoing introgression, we cannot rule out the possibility that more complex, untested models might provide a better fit to the data (see 'Methods' section). This inferred strict genetic isolation despite close relatedness and reproductive interference between M. achilles and M. helenor nevertheless suggests that ecological factors might prevent gene flow in sympatry.
Temporal segregation between sympatric sister species. The analysis of the temporal variations of flight activity in a markrecapture experiment revealed a striking difference in patrolling time among species (Kruskal-Wallis test: Chi-square = 179.7, P < 0.001, df = 3), with little overlap between the sister species M. achilles and M. helenor (Fig. 4). Males M. helenor patrolled earlier than M. achilles (mean patrolling time ± s.d. = 11:14 ± 00:45 vs. 12:35 ± 00:40, respectively). Patrolling time in M. deidamia (12:40 ± 00:46) however overlapped with M. achilles. Time of (re) capture of the same individual in M. achilles was notably stable (correlation between time of first vs. second capture: r = 0.40; P = 0.05), suggesting a regularity in patrolling time at the individual level in this species (Fig. 4). Whether such individual temporal regularity is genetically determined or reflects a plastic behaviour 34,35 remains to be investigated. Reproductive and/or aggressive interference are generally expected to promote spatial or temporal habitat segregation between species because it reduces the cost of negative interspecific interactions 28,36 . The observed shift in flight activity peaks suggests that reproductive interferences among the morphologically convergent Morpho species have favoured the evolution of divergent temporal niches in M. achilles and M. helenor. Remarkably, this temporal segregation extended to the more distantly-related Morpho menelaus, that flies earlier in the morning (Fig. 4). In particular, no overlap between its activity range and that of M. achilles was detected, mitigating the potentially strong interference between these species suggested by the attraction of patrolling M. achilles towards M. menelaus dummy (Fig. 1b). Whether the temporal segregation observed between these sympatric Morpho species was involved in a sympatric speciation process or stem from reinforcement after secondary contact between these closely-related species remains to be elucidated, as well as the genetic vs. plastic origin of this behavioural trait.
Temporal segregation in response to resource competition or reproductive interference has been reported in other taxa [37][38][39] , but mostly at a larger temporal scale (e.g. seasonal segregation). Peaks of patrolling activity between sister Morpho species are shifted by less than 2 h, representing a remarkably fine temporal scale for a shift in reproductive behaviour 39 . Nonetheless, comparable fine-scale temporal partitioning of sexual activity has been observed in other moth and butterfly species 35,40 , suggesting that it may be relatively common in Lepidoptera.
Altogether, our results suggest that there is a potential for strong interspecific competition among males from Morpho species in sympatry, but that this competition is mitigated at least in part by their temporal segregation. Heterospecific courtships may also be limited, if the temporal segregation observed in  females' activities. This could not be tested here as females were almost never seen flying along the riverbank: future studies specifically exploring females behaviour should thus be performed. Synchronisation of mating activity between sexes is however likely 41,42 , because any deviation of males relative to females activity would reduce the probability of intraspecific mating 34,35 . In contrast, a shift in mating time among species would act as a powerful isolation mechanism 39 , explaining how co-mimetic Morpho species can maintain a similar colour pattern in the same habitat while remaining sexually isolated. Divergence in daily phenology in butterfly mating activities may be a widespread process enabling the persistence of diversityrich assemblages, as suggested by several reports of temporally structured sexual activities in other butterflies [43][44][45] , including closely-related species 46 . Our study shows that the coexistence of closely-related species can generate complex ecological interactions, both mutualistic (mimicry) and antagonistic (reproductive interference), that could be mitigated by shifts in temporal niches. The evolution of multiple traits, morphological and behavioural, may thus favour species diversification in sympatry by partitioning niche along different ecological dimensions.

Methods
Study site and population. The study was conducted between July and October 2019 in the North of Peru. We focused on populations of coexisting Morpho species present in the regional park of the Cordillera Escalera (San Martin Department) near the city of Tarapoto. Both the capture-recapture and the dummy experiment were performed at the exact same location, on the bank of the Shilcayo river (06°27′14.364″S, 76°20′45.852″W).
DNA extraction and RAD-Sequencing. Thirty-one wild males caught on the study site were sequenced to perform population genomic analyses (M. achillesn = 13, M. helenor-n = 10 and M. deidamia-n = 8). DNA was extracted from each sample from a slice of the thorax, using Qiagen kit DNeasy Blood & Tissue. DNA quantification (using the microfluorimetric method) and quality controls (using electrophoresis and spectrophotometric method) were performed prior to sequencing. RAD-library preparation and sequencing were performed at the MGX-Montpellier GenomiX platform (Montpellier, France). DNA was digested with the Pst1 enzyme and the library was prepared according to Baird and Etter' protocol 47 in a slightly modified version. Paired-end RAD-sequencing was performed on a 2 lanes flow cell of an Illumina HiSeq2500 in a rapid mode so that reads (125 bp) were expected to be of high quality with no missing base (N content). We obtained 299 million sequences, comprising R1 and R2 reads for each sequenced fragment. Adapters were removed from the reads.
Read quality control, alignment and dataset generation. Read quality was assessed with FastQC v0.11.9 (http://www.bioinformatics.babraham.ac.uk/projects/ fastqc/). The per base sequence quality was high across all reads (no lower than 36 for R1 and 32 for R2) with an average quality score of 39 (40 being the maximum). Overall, FastQC highlighted the high quality of the sequencing data, allowing us to skip the step of read trimming.
The data were demultiplexed, assigning each sequence to its sample ID and the reads were aligned using Stacks V2.5 (http://catchenlab.life.illinois.edu/stacks/) 48,49 . Parameters were set following the 80% polymorphic (r80) loci rule, which only considers loci shared by at least 80% of the samples 50 . The optimised parameters are 'max distance between stacks' (inside each sample) and 'number of mismatches between stacks' (between samples). Every other parameter was kept to default values. After aligning all reads, we selected 2740 biallelic loci shared by all samples, including 88,889 SNPs in total. Each locus had a length of 463.12 bp on average (range [343; 908]). These loci are assumed to be evenly distributed throughout the genome but cover only a limited portion of the genome (around 0.5%). Datasets were stored in a VCF file (containing all the SNPs found in the alignment) and a fasta file (containing the two alleles found at every locus for each sample). To run DILS-ABC inferences, Stacks fasta files were converted to another fasta format compatible with DILS (https://github.com/CoBiG2/RAD_Tools).
Demographic inferences. Eight categories of demographic models were compared, according to temporal patterns of introgression. This was done to answer two questions on gene flow in Morpho: (1) is there ongoing migration between M. helenor and M. achilles? (2) do M. helenor and/or M. achilles exchange alleles with M. deidamia? This was assessed by an ABC approach using a version of DILS adapted to samples of three populations/species 32 . Since Stacks does not report monomorphic RAD loci, the ABC analysis was conditioned in the same way, by excluding monomorphic loci from the simulations. Focusing on polymorphic loci may only limit our ability to estimate the absolute values of parameters (i.e. population sizes expressed in numbers of individuals, and ages of past events expressed in numbers of generations); nevertheless, this framework excluding monomorphic loci still allows reliable comparisons of models 51 and estimations of relative parameter values, as performed to investigate the human history 51 .
A generalist model was studied (Supplementary Fig. 12). This model describes an ancestral population subdivided in two populations: the ancestor of M. deidamia and the common ancestor of M. helenor/M. achilles. The latter population was further subdivided into the three species/populations currently sampled. Each split event is accompanied by a change in demographic size, the value of which is independent of the ancestral size. In addition, given clear genomic signatures for recent demographic changes with largely negative Tajima's D, we implemented variations for the effective sizes of the three modern lineages at independent times. Finally, migration can occur between each pair of species/ populations. Migration affecting the M. helenor/M. achilles pair can either be the result of secondary contact after a period of isolation (ongoing migration), or of ancestral migration (current isolation) as in 50,52 .
As this model is over-parameterised, our general strategy is to investigate the above two questions by comparing variations of this generalist model. Thus, to test the gene flow between M. helenor and M. achilles, we compared two categories of models. (1) With random parameter values for all model parameters including the ongoing migration between M. helenor and M. achilles (gene flow resulting from a secondary contact between them); (2) as above, but with the migration between M. helenor and M. achilles set to zero after a randomly drawn number of generations following their split. An overlap between 'current isolation' and 'ongoing migration' models can occur when the transition time (from ancestral migration to current isolation forward in time for a 'current isolation' model; or from ancestral isolation to ongoing migration forward in time for an 'ongoing migration' model) tends towards the extreme values 0 or T split hel-ach ( Supplementary Fig. 12). To reduce this effect, the transition times were drawn in a Beta distribution with parameters (α = 5, β = 1) when migration has to be restricted to a past period, and in a Beta distribution with parameters (α = 1, β = 5) when migration is assumed to occur after a recent secondary contact.
When two broad categories of models are statistically compared, each category is represented by simulations performed under the four sub-models allowing or not allowing genomic heterogeneities for effective sizes (Ne) and for migration rates (N.m). For instance, to test for gene flow between M. helenor and M. achilles, the model of 'ongoing migration' is actually represented by simulations with the four possible combinations of homogeneity/heterogeneity, all labelled as being 'ongoing migration'.
As for any inferential analysis, it is important to recognise that the bestsupported model is based on a classification of models within a studied set. Intermediate models, with more subtle cycles of genetic isolation and secondary contact could produce a better fit to the data, but it would be surprising to detect a strong support for the model assuming a lack of recent gene flow, if the most recent secondary contact of such cyclicity induced elevated gene flow.
For each model, 50,000 simulations using random combinations of parameters were performed. Parameters were drawn from uniform prior distributions. Population sizes were sampled from the uniform prior [0-1,000,000] (in diploid individuals); the older time of split was sampled from the uniform prior [0-8,000,000] (generations); ages of the subsequent demographic events were sampled in a uniform prior between 0 and the sampled time of split. Migration rates 4.N.m were sampled from the uniform prior [0-50]. Both migration rates and effective population sizes are allowed to vary throughout the genomes as a result of linked selection, following refs. [53][54][55] .
On each simulated dataset, we calculated a vector of means and standard deviations for different summary statistics: intraspecific statistics (π for M. helenor, π for M. achilles, π for M. deidamia, θ W for M. helenor, θ W for M. achilles, θ W for M. deidamia, Tajima's D for M. helenor, Tajima's D for M. achilles, Tajima's D for M. deidamia) and interspecific statistics (gross divergence, net divergence and F ST for all three possible pairs; ABBA-BABA D). Our version of DILS includes part of the DaDi 56 and Moments 57 strategy involving the identification of the best model proposed demographic model from the molecular patterns of polymorphism and divergence (proportion of shared polymorphisms, fixed differences between species, exclusive polymorphisms, etc.), excluding monomorphic loci. Thus, only loci containing at least one SNP in an alignment of the three species studied are considered, including singletons. Importantly, each locus carrying at least one SNP in a tri-specific alignment is associated with a mutation rate assumed to be 3 · 10 −9 mutations per generation and per base pair to convert demographic parameters into demographic units from coalescence units.
We first conditioned the mutations occurring during coalescent simulations by using theta (=4 · N · µ · L i ; where N is the effective population size, µ the mutation rate per nucleotide and per generation; L i the length of locus i). The number of simulated segregating sites for a given locus strongly depends on the coalescent history (i.e the total length of the simulated coalescent tree), occasionally generating monomorphic loci. To confirm that the inferences are not impacted by differences in the number of monomorphic loci in the simulated datasets, we then used an alternative simulation approach, by randomly placing in simulated coalescent trees a fixed number of mutations corresponding to the observed number of SNPs for each locus. Thus, a randomly simulated dataset consists of 2740 loci whose lengths (ranging from 339 to 894 nucleotides) and number of SNPs (ranging from 1 to 91) individually match the properties of the observed loci in the actual dataset. Since the results drawn from both approaches were similar, we report only the estimations provided by the simulations based on the actual number of SNPs. Comparisons between the two approaches can be found in supplementary (Supplementary Tables 8,9).
Statistical comparisons between simulated and observed statistics were performed using the R package abcrf version 1.8.1 58,59 .
Mark-recapture experiment. To estimate the timing of patrolling activity among Morpho species, we performed capture-mark-recapture between 9 a.m. and 2 p.m. (flight activity in Morpho is drastically reduced in the afternoons at this site) during 17 sunny days. Although on a few days, capture was cancelled because of bad weather annihilating butterfly activity, the 17 capture sessions were mostly consecutives, as they were performed in a 22 days period (Supplementary Table 1 and Supplementary Fig. 15). All butterflies were captured with hand-nets, identified at the species level, and numbered on their dorsal wing surface using a black marker. The exact time of each capture was annotated. Butterflies captured while inactive, such as those laying on a branch or on the ground were excluded from the analysis to focus exclusively on actively patrolling individuals. We measured patrolling time for a total of 295 occasions, including 78 recaptures (i.e. 217 individuals were captured at least once). All captured individuals were males. Individuals M. achilles were the most frequently captured (n = 121), followed by M. helenor (n = 95). Individuals M. deidamia were about half less captured (n = 48), and individual M. menelaus were the least captured (n = 34). Because striking differences in patrolling time were observed among Morpho species, we used time of the day as a predictor of species identity in order to distinguish between M. helenor and M. achilles in the below-described experiment because butterflies from these two species are morphologically too similar to be identified while flying ( Supplementary  Fig. 13). After the 17 nearly-consecutive days of capture, one day of capture was repeated every 2 weeks during 2 months in parallel to the dummy experiment (described below), to verify that temporal activity was stable over time (Supplementary Fig. 13).
Estimating population size from mark-recapture data. Based on capturerecapture histories, we estimated individual abundance for each species using a loglinear model implemented in the R package Rcapture version 1.4.3 60 (Supplementary Fig. 15). Given the short duration the sampling period (22 days) relative to the longevity of adult Morpho butterflies (several months 61 ), we used a closedpopulation model assuming no effect of births, deaths, immigration and emigration. Abundance was estimated in Morpho helenor and M. achilles only, as capture and recapture events were too few in the other species (M. deidamia and M. menelaus) to allow estimating population size (Supplementary Table 1).
Experiment with dummy butterflies. We investigated the response of patrolling males to sympatric conspecifics, congeners and of exotic conspecifics, using dummies placed on their flight path. Dummies were built with real wings dissected and washed with hexane to remove volatile compounds and cuticular hydrocarbons, ensuring to test only the visual aspect of the dummies. We mounted the wings on a solar-powered fluttering device (Butterfly Solar Héliobil R029br) that mimics a flying butterfly, thereby increasing the attractiveness of the dummy. The fluttering dummy was positioned on the riverbank, and placed at the centre of a 1 m 3 space delimitated with four vertical stacks (Fig. 1a). The set-up was continuously monitored by a human observer and filmed using a camera (Gopro Hero5 Black set at 120 images per second) mounted on a tripod. Patrolling Morpho butterflies that deviated from their flight path to approach the dummy but did not enter the cubic space were categorised as approaching. Any Morpho butterfly entering the cubic space was considered as interacting with the dummy. Those passing without showing interest to the setup were categorised as passing. The category of behaviour and the exact time of the butterfly responses were annotated on site by the human observer. Patrolling individuals were mainly identified at the species level by the observer on the site: M. menelaus can be easily distinguished from M. deidamia, and these two species are also quite different from M. helenor and M. achilles. However, the sister species M. helenor and M. achilles cannot be discriminated during flight, and we thus rely on an indirect method, based on flight hours, to infer the species identity of wild visitors looking as a M. helenor/M. achilles (Supplementary Fig. 13). Note that removing data with the highest levels of uncertainty in species identity (i.e. when discarding visits performed in the period where M. helenor/M. achilles temporally overlap) does not quantitively affect our results ( Supplementary Fig. 14 and Supplementary Tables 5, 6). Using the recorded video, we also measured the duration of the interactions (i.e. the time spent in the cubic space) occurring between patrolling male and the dummy. The ten dummies were each tested during 4 sunny days from 9 a.m. to 2 p.m. (i.e. during 5 h). This resulted in 40 days of experiment over which each dummy was left fluttering on the river bank for a combined duration of 20 h. Dummies were randomly attributed to each day of the experiment. Mark-recapture data suggested a very low rate of individuals passing through the site several times per day (mean percentage of recapture within the same day = 0.95%), thus limiting potential pseudoreplication within each dummy replicate. We recently showed that intraspecific variation in wing colour pattern within the locality is very low in these species 25 . Using a single dummy per sex and species, as done here, should thus have little impact on the observed behaviours.
In order to control for variation in weather (affecting both the activity of patrolling butterflies and of the solar-powered device), we collected hourly data on the percentage of cloud cover for the period and location of our experiment (available at https://www.visualcrossing.com/). A percentage of cloud cover was then associated with all the behavioural observations, and used as a control variable in all statistical analyses.
Three-dimensional kinematics of flight interaction with the dummies. To test whether Morpho males showed different flight behaviours when interacting with the male and female dummy, we filmed the flight interactions using two orthogonally positioned video cameras (Gopro Hero5 Black, recording at 120 images per second) around the dummy setup (Fig. 1a). Stereoscopic video sequences obtained from the two cameras were synchronised with respect to a reference frame (here using a clapperboard). Prior to each filming session, the camera system was calibrated with the direct linear transformation (DLT) technique 62 by digitising the positions of a wand moved around the dummy. Wand tracking was done using DLTdv8 63 , and computation of the DLT coefficients was performed using easyWand 64 . After spatial and temporal calibration, we also used DLTdv8 to digitise the three-dimensional positions of both the visiting (real) butterfly and the dummy butterfly at each video frame by manually tracking the body centroid in each camera view. Butterfly positions throughout the flight trajectory were postprocessed using a linear Kalman filter 65 , providing smoothed temporal dynamics of spatial position, velocity and acceleration of the body centroid. Based on these data, we investigated how spatial position, speed and acceleration of the visitor butterfly varied over the course of the interaction. We proceeded by dividing space into 10 cm spherical intervals around the dummy position ranging from 0 to 1.2 m distance (this step standardises interactions of different durations), and computed the proportion of time spent, the mean speed and acceleration of the interacting butterfly within each distance interval (Fig. 2). We analysed a total of 28 interactions performed by individual Morpho achilles male, including 14 with the dummy of its conspecific male and 14 with the dummy of its conspecific female. Analysed interactions lasted in average 1.44 ± 0.87 (mean ± sd) s.
Statistical analysis of behavioural experiments. Differences in patrolling time were assessed by testing the effect of species on time of capture using Kruskal-Wallis test. To test the effect of visitor identity and dummy characteristics on the number of approaches and interactions, we performed logistic regressions. Approach was treated as a binary variable, where 0 meant 'passing without approaching' and 1 meant 'approaching the dummy setup'. For the interactions, we only considered individuals approaching the setup, such as 0 meant 'approaching without entering the cubic space' and 1 meant 'entering the cubic space'. This allowed getting rid of the uncertainties on whether passing individuals had actually seen the setup or not. We first tested the effect of visiting species on approach and interaction while controlling for dummy's characteristics to test for intrinsic differences in territoriality (or 'curiosity') among species. We then tested the effect of the dummy sex and identity on approach and interaction separately in Morpho helenor and M. achilles. The percentage of cloud cover was also included in the models to control for variation in dummy movements (generated by the solarpowered device), potentially affecting the butterfly response (Supplementary  Tables 3 and 4). We further tested if variation in wing area and proportion of iridescent blue among dummies affected the frequency of approach and interaction, again using logistic regression analyses ( Supplementary Fig. 7). Statistical significance of each variables was assessed using likelihood ratio tests comparing logistic regression models 66 . Finally, we tested the effect of dummy sex and identity on the duration of interaction using Kruskal-Wallis tests.
Based on the flight kinematic data, we investigated whether flight behaviour during the interaction differed with male vs. female dummies. We ran a mixedeffects model testing the effect of (1) the sex of the dummy and of (2) the distance from dummy (fixed effects), on the proportion of time spent (fixed effects), using the flight ID as a random effect. The flight ID corresponds to the behaviour of a single wild males flying within the 'interaction space'. Specifically, we tested for the statistical interaction between the sex of the dummy and distance from dummy on the proportion of time spent in the different distance intervals. We then similarly tested for difference in acceleration over the course of the flight interaction, by testing the effect of (1) the sex of the dummy and of (2) the distance from dummy (fixed effects), on the acceleration, with the flight ID as a random effect. We focused on the statistical interaction between the sex of the dummy and the distance from dummy on the mean acceleration in the different distance intervals.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
All data supporting the findings of this study are provided as a Source Data file. The RNA-seq data generated in this study have been deposited to the National Center for Biotechnology Information under accession number PRJNA739839 (https:// www.ncbi.nlm.nih.gov/). Source data are provided with this paper.

Code availability
Genomic analyses were performed using the statistical analysis platform DILS (https:// eep.univ-lille.fr/en/productions-2/dils-software/). A modified version of DILS adapted for studying divergence among three species/populations was used for this study. Source codes can be freely used from GitHub: https://github.