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. achilles—n = 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’ protocol47 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 samples50. 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/species32. 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 models51 and estimations of relative parameter values, as performed to investigate the human history51.
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 in50,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 Tsplit 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 best-supported 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 FST for all three possible pairs; ABBA-BABA D). Our version of DILS includes part of the DaDi56 and Moments57 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 · µ · Li; where N is the effective population size, µ the mutation rate per nucleotide and per generation; Li 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.158,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 capture-recapture histories, we estimated individual abundance for each species using a loglinear model implemented in the R package Rcapture version 1.4.360 (Supplementary Fig. 15). Given the short duration the sampling period (22 days) relative to the longevity of adult Morpho butterflies (several months61), we used a closed-population 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 m3 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 species25. 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) technique62 by digitising the positions of a wand moved around the dummy. Wand tracking was done using DLTdv863, and computation of the DLT coefficients was performed using easyWand64. 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 post-processed using a linear Kalman filter65, 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 solar-powered 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 models66. 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 mixed-effects 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.
Source: Ecology - nature.com
