Cryptic speciation associated with geographic and ecological divergence in two Amazonian Heliconius butterflies

The evolution of reproductive isolation via a switch in mimetic wing coloration has become the paradigm for speciation in aposematic Heliconius butterflies. Here, we provide a counterexample to this, by documenting two cryptic species within the taxon formerly considered Heliconius demeter Staudinger, 1897. Amplified fragment length poly-morphisms identify two sympatric genotypic clusters in northern Peru, corresponding to subspecies Heliconius demeter ucayalensis H. Holzinger & R. Holzinger, 1975 and Heliconius demeter joroni ssp. nov. These subspecies are reciprocally monophyletic for the mitochondrial genes COI and COII and the nuclear gene Ef1 α , and exhibit marked differences in larval morphology and host plant use. COI sequences from 13 of the 15 currently recognized subspecies show that mtDNA differences are reflected across the range of H. demeter , with a deep phylogenetic split between the southern and northern Amazonian races. As such, our data suggest vicariant speciation driven by disruptive selection for larval performance on different host plants. We raise Heliconius demeter eratosignis (Joicey & Talbot, 1925) to Heliconius eratosignis based on nomenclatural priority, a species also comprising H. eratosignis ucayalensis comb. nov. and three other southern Amazonian races. Heliconius demeter joroni spp. nov. remains within H. demeter s.s. , along with northern Amazonian and Guianan subspecies.


INTRODUCTION
Cryptic species can be defined as species that are, or have been, erroneously classified as a single nominal species due to their superficial morphological similarity (Bickford et al., 2007). In recent years, the integration of DNA sequences into taxonomy has led to the recognition of increasing numbers of such species (Hebert et al., 2004). Whether the origins of cryptic species can be consistently ascribed to any particular evolutionary processes is unclear (Bickford et al., 2007). On the face of it, this might seem unlikely, given that they are defined principally through humans' visual perception of morphology. For example, while divergent selection may often lead to sister species with markedly different body shapes or colours (Jiggins et al., 2001b;Langerhans et al., 2007), in other cases it may affect traits causing reproductive isolation that have no clear morphological basis, such as behaviours (Janzen et al., 2009). Nonetheless, some theories of speciation might be more predisposed to the creation of cryptic species, such as when reproductive isolation results from the chance fixation of different, epistatic incompatibilities in separate populations subject to similar selective pressures (Clarke et al., 1988;Mani & Clarke, 1990;Orr, 1995;Turelli & Orr, 2000). It also has been shown that, at least for butterflies in the western Mediterranean, cryptic species are rarely sympatric (Vodă et al., 2015). While this may reflect to some extent the ability of taxonomists to diagnose species with limited range overlap, it is also consistent with phenotypic similarity constraining coexistence (Pigot & Tobias, 2013).
Heliconius butterflies are chemically defended and aposematic, i.e. they advertise their defence to wouldbe predators using bright colours on their wings. To minimize the per capita cost incurred while predators learn the association between the warning signal and prey unprofitability, many Heliconius species mimic one another. This mutualistic interaction is known as Müllerian mimicry (Müller, 1879). Within Heliconius, a number of distinct mimetic phenotypes exist (e.g. blue and yellow, red and black patterns). Groups of sympatric species exhibiting the same phenotype are said to be co-mimics in a 'mimicry ring'. It has been convincingly shown that when a population switches to a different mimicry ring, this can contribute to reproductive isolation. The reasons for this are twofold. First, hybrids with intermediate colour patterns are selected against by predators that do not recognize them as aposematic (Merrill et al., 2012). Second, Heliconius males have been shown to preferentially court females with similar colour patterns to their own (Jiggins et al., 2001b;Jiggins et al., 2004;Kronforst et al., 2006;Mavárez et al., 2006;Chamberlain et al., 2009;Merrill et al., 2011), even when those females belong to a different species (Estrada & Jiggins, 2008). Because divergence in an ecologically relevant adaptive trait also creates reproductive isolation, Heliconius have become a prime example of so-called 'ecological speciation' (Nosil, 2012).
Nonetheless, there are instances of Heliconius sister species that do not appear to have diverged in wing colour pattern. For example, Heliconius sara (Fabricius, 1793) and H. leucadia Bates, 1862 are sympatric sister species with almost identical blue and yellow phenotypes. Heliconius numata (Cramer, 1780) and H. ismenius Latreille, 1817 are parapatric sister species with similar 'tiger' colour patterns. In addition, modern taxonomy and DNA sequencing have revealed a number of cryptic races belonging to the H. cydno/timareta superspecies from the tropical eastern Andes (Brower, 1996;Lamas, 1997;Giraldo et al., 2008;Mallet, 2009;Mérot et al., 2013;Arias et al., 2017). These taxa were hitherto unrecognized as members of the H. cydnotimareta clade because they exhibit colour patterns extremely similar to those of sympatric subspecies of H. melpomene (Linnaeus, 1758), itself the sister to the H. cydno-timareta lineage. In some cases, this striking phenotypic similarity is likely due to adaptive introgression of colour patterns between H. melpomene and H. timareta Hewitson, 1867(Heliconius Genome Consortium, 2012. These examples suggest that speciation in Heliconius may sometimes occur without a mimicry shift, and demonstrate that closely related, co-mimics can maintain their identities in sympatry, despite occasional hybridization (Mérot et al., 2017).
Based on previous systematic research, H. demeter Staudinger, 1897 was held to comprise 15 described subspecies with red, yellow and black phenotypes ( Fig. 1) that participate in the 'dennis-rayed' Heliconius mimicry ring (Brown & Benson, 1975;Lamas, 2004). The taxon is widely distributed throughout most of Amazonia and the Guiana shield, but is usually scarce when compared to closely related co-mimics, such as H. erato (Linnaeus, 1758). Interestingly, several north Amazonian and Guianese races of H. demeter are sexually dimorphic, with hindwing rays in males fused at their base to form a bar. Sexual dimorphism in colour pattern is rare in Heliconius, and only one other species exhibits it prominently: Heliconius nattereri C. Felder & R. Felder, 1865, from south-eastern Brazil.
In the easternmost cordillera of the Andes in northern Peru, we discovered what we at first took to be two H. demeter races, H. demeter cf. demeter and H. demeter cf. ucayalensis H. Holzinger & R. Holzinger, 1975, flying together near the city of Tarapoto. Heliconius demeter cf. demeter is sexually dimorphic, but H. demeter cf. ucayalensis is not. At first we viewed these taxa as somewhat divergent subspecies, since there are contact zones between many butterfly subspecies in this area (Dasmahapatra et al., 2010). In the Downloaded from https://academic.oup.com/zoolinnean/advance-article-abstract/doi/10.1093/zoolinnean/zly046/5066664 by University of York user on 06 August 2018 present study, we show that these sympatric subspecies of 'H. demeter' in fact comprise two distinct species, corresponding to H. demeter cf. demeter and other northern and central Amazonian subspecies, and H. demeter cf. ucayalensis and the south Amazonian races. In accordance with nomenclatural priority, the southern clade is recognized as H. eratosignis (Joicey & Talbot, 1925), a species comprising four subspecies (Lamas & Jiggins, 2017;Supporting Information, Table S1), and this nomenclature is adopted in this paper from here on. Additionally, it was noted that H. demeter cf. demeter specimens from Tarapoto are divergent from those in the H. demeter type locality near Iquitos, and accordingly this population is here described as a new subspecies: Heliconius demeter joroni ssp. nov.

MORPHOLOGICAL AND BEHAVIOURAL ANALYSIS
To identify species-specific diagnostic characters in the 15 currently recognized subspecies of H. demeter and H. eratosignis, all type series and specimens held in the Natural History Museum London (NHMUK) were examined. In addition, we examined holotypes, allotypes, syntypes and other material held at the Florida Museum of Natural History (FLMNH), the Museum für Naturkunde, Berlin (MNB), the Natural History Museum at the San Marcos National University, Lima, Peru (MUSM), the Naturhistorisches Museum, Wien (NHMW), the National Museum of Brazil, Rio de Janeiro (MNRJ), the Museum of Zoology 'Adão José Cardoso' at the University of Campinas, Brazil (ZUEC) and the Museum of Zoology at the University of São Paulo, São Paulo, Brazil (MZUSP) (Supporting Information, Table S2).
For morphometric analyses of wing shape, images of the ventral and dorsal surfaces of dissected forewings and hindwings of 75 H. eratosignis ucayalensis, 31 H. demeter joroni ssp. nov. and 16 H. demeter bouqueti Nöldner, 1901 specimens from Tarapoto and French Guiana were captured using either a high-resolution flatbed scanner or a Nikon D90 digital camera with a Nikon micro 105/2.8GEDVR lens. In addition, we conducted a global geometric morphometric analysis using 31 photographs of museum specimens representing eight other subspecies. All specimens used in morphometric analysis are shown in Supporting Information,  Turner, 1966 and H. d. terrasanta appear to conform to the type specimens only around the type localities (in Terrasanta, Pará, and in Guyana). Between these, most populations appear to be either polymorphic or exhibit intermediate phenotypes (mixed square and cross symbols in the map). Heliconius demeter ssp. nov. refers to three males in the FLMNH recognized by W. Neukirchen as distinct from other described subspecies. These individuals may prove to have affinities to H. demeter titan.
Forewing and hindwing shape were described using 20 and 18 landmarks, respectively, which were placed at vein intersections and vein termini on the ventral side (Supporting Information, Fig. S1). Standard tests of repeatability were carried out by taking the landmarks five times per wing on subsamples of five butterflies from a single subspecies and sex. Landmark coordinates were digitized using TpsDig2 (Rohlf, 2010) and superimposed using a general Procrustes analysis (Bookstein, 1991;Zelditch et al., 2004). Wing size was measured using the log-transformed centroid size (Bookstein, 1991). Differences in size between H. d. joroni ssp. nov. and H. e. ucayalensis were investigated with a one-way ANOVA, with size as the response, and species and sex as predictive factors. P-values were corrected for multiple comparisons following Benjamini & Hochberg (1995).
To study shape, dimensionality reduction was employed to correct for the effect of using a large number of variables relative to the number of specimens. We used the minimum subset of principal components (PCs) that minimized the total cross-validated misclassification percentages between groups defined a priori (Baylac & Friess, 2005). To explore shape differences between H. d. joroni ssp. nov. and H. e. ucayalensis, a MANOVA was applied to the PC subsets, with shape as the response and sex and species as predictive factors. Given the high sexual dimorphism, species discrimination based on shape was investigated for each sex separately through a Canonical Variate Analysis (CVA), with a leave-one-out cross-validation procedure (CV). All statistics and morphometrics were performed in R 2.13.1 (R Development Core Team, 2011) with ade4 (Chessel et al., 2004) and Rmorph libraries (Baylac, 2007).
Genitalia of three male H. d. joroni ssp. nov. and seven male H. e. ucayalensis collected from Tarapoto were prepared from material preserved in salt-saturated DMSO. The tips of the abdomens were removed and soaked in 10% KOH for 10 min at 70 °C, and then transferred to distilled water. The scales were first removed with a fine brush and the valves extruded. The genitalia were then removed and further cleaned. Temporary slides were prepared in 25% ethanol, and the interior surfaces of each left valva were photographed.
Observations on host plant use and larval morphology were made near Tarapoto, Peru. To supplement field observations of host plant use, wild caught adult females were placed in a cage with 22 locally common Passiflora species (Supporting Information, Table S4) and allowed to oviposit. Geographic localities for H. demeter and H. eratosignis were obtained from those published in (Rosser et al., 2012) and supplemented with subsequent collections by NR in Bolivia, Brazil, French Guiana, Peru and Suriname between 2011, by AVLF in Mato Grosso and Acre from 1994to 2016, and by Keith Brown (from 1970to 1999) and Eurides Furtado (from 1978to 1998 in Brazil.

MOLECULAR ANALYSIS
Details of the specimens used for molecular work are shown in Supporting Information, Table S5. Wings were removed from samples collected in French Guiana and around Tarapoto, and the bodies preserved at -20 °C in salt-saturated DMSO. Both wings and tissue of the French Guiana specimens are held at the Muséum National d'Histoire Naturelle, Paris (MNHN), while the Peruvian specimens are held at the University of York, UK. In addition, single legs representing 11 other subspecies were obtained from dried museum specimens in the FMNH (identification numbers beginning with 'KW' in Supporting Information, Table S5). DNA was extracted from these legs using the QIAamp DNA Micro Kit (QIAGEN), and from onethird of the thorax of the remaining specimens using the DNeasy Blood and Tissue Kit (QIAGEN).
Approximately 2200 bp of mtDNA comprising cytochrome oxidase I (COI), tRNA-leu and the 5' end of cytochrome oxidase II (COII) was amplified by PCR in three sections, for seven H. demeter joroni ssp. nov., 12 H. eratosignis ucayalensis and 12 H. demeter bouqueti. Four autosomal nuclear genes Elongation factor 1-α (Ef1α), Tektin, Ribosomal protein L5 (Rpl5), Mannose-phosphate isomerase (Mpi) and the sexlinked Triose phosphate isomerase (Tpi) were also successfully sequenced for varying numbers of these three taxa. Only small amounts of degraded DNA could be obtained from the museum specimens. Therefore, for these samples the first ~760 bp of COI was amplified in two shorter sections. All PCR products were cleaned and cycle-sequenced with the PCR primers using the BIG DYE TERMINATOR v.3.1 Cycle Sequencing Kit (Applied Biosystems), and sequences obtained using an ABI3730xl DNA Analyzer (Applied Biosystems). Supporting Information, Table S6 contains details of the primers used and PCR conditions. Sequences from Heliconius species in the sara-sapho clade were downloaded from Genbank to act as outgroups. GenBank accession numbers for all sequences used are provided in Supporting Information, Table S7 and Table S8. Sequences were aligned using ClustalW, and the alignments then checked by eye.
Phylogenetic analysis of sequence data was carried out using MEGA7 (Kumar et al., 2016). For each gene, we found the nucleotide substitution model that best described the substitution pattern using all sites, a Neighbour-joining tree and Bayesian Inference Criterion (BIC). We then found the maximum likelihood (ML) tree for each gene assuming the selected model of sequence evolution, and estimated node Downloaded from https://academic.oup.com/zoolinnean/advance-article-abstract/doi/10.1093/zoolinnean/zly046/5066664 by University of York user on 06 August 2018 support using 1000 bootstrap replicates. To obtain higher resolution nuclear data, eight specimens each of H. eratosignis ucayalensis, H. demeter joroni ssp. nov. and H. demeter bouqueti were genotyped using four AFLP primer combinations: TaqI-CAG with EcoRI-ATG, TaqI-CGA with EcoRI-AGC, TaqI-CAG with EcoRI-AGC and TaqI-CCA with EcoRI-ACA. The AFLP protocol used is similar to that described in Vos et al. (1995), and the primer sequences and reaction conditions are described in Madden et al. (2004). The AFLP products were resolved by electrophoresis through 6% acrylamide gels, visualized by autoradiography, and scored by eye. A total of 81 loci were polymorphic and could be scored unambiguously.
The Bayesian clustering program STRUCTURE 2.2 (Pritchard et al., 2000) was used to evaluate the number of genetic clusters indicated by these AFLP genotypes, using standardized inference criteria (Evanno et al., 2005). Following a 100 000-step burnin period, data were collected over 100 000 Markov chain Monte Carlo repetitions. STRUCTURE analysis was carried out on the dataset, increasing K from 1 to 10. At each value of K, the analysis was repeated three times to check between-run consistency. The AFLP data were also used to calculate pairwise Nei-Li genetic distances (Nei & Li, 1979) between all genotyped individuals. These distances were then used to calculate average genetic distances between each of the three taxa.

ADULT MORPHOLOGY
Examination of H. demeter joroni ssp. nov. and H. eratosignis ucayalensis wings from Peru revealed two diagnostic morphological characters strictly concordant with the mitochondrial DNA classification of these two species.
(1) In the proximal region of the narrow costasubcosta space on the underside of the forewing, H. demeter joroni ssp. nov. exhibits a strong yellow 3-5mm-long streak placed in the anterior half of the space along the costal vein, often associated with black scales posteriorly ( Fig. 2A). In H. eratosignis ucayalensis this region is uniformly orange (Fig. 2B). Brown & Benson (1975) also noticed this character difference between northern and southern Amazonian populations, but did not recognize its significance, probably because they lacked a long series of these taxa from a sympatric population.
(2) In males of H. demeter joroni ssp. nov., the red rays on the dorsal hindwing fuse to form a hindwing bar (Fig. 2C), while in males of H. e. ucayalensis they do not (Fig. 2D). This character is inapplicable to females, all of which have unfused red rays, and to the geographic forms of H. demeter from north-eastern South America that lack rays. One clear difference in genital morphology was observed between the males of the H.  this region has a characteristic convex depression Supporting Information, Fig. S2). However, the utility of this trait is unclear given the small sample sizes.
Using the presence/absence of the yellow costal streak on the ventral forewing, the existing 15 named subspecies of H. demeter could be unambiguously classified as either belonging to H. demeter or H. eratosignis COI haplogroups (see below). With the exception of Heliconius demeter titan Neukirchen, 1995, in all male specimens the presence of the yellow costal streak was also perfectly concordant with fused or reduced hind wing rays (Supporting Information, Hindwing: On dorsum the grey friction patch is narrow, and the ray in cell Rs-M 1 is strongly present, forming the anterior tip of the bar of fused rays (in H. e. ucayalensis rays are unfused and the friction patch is broad, leading to almost complete loss or reduction to a smudge of the ray). On the ventral side a yellow costal streak, a single row of white submarginal dots along the anal margin and some diffuse red spots at the bases of Cu 2 , Sc+R 1 and the discal cell. Rays reduced relative to the dorsal side, and unfused.

Female
Forewing: Length: 35-39.5 mm, mean = 36.8 mm, N = 5. As the male, except no friction patch or greenish tinge to forewing postmedian band on dorsum, and no greenish spot in the middle of cell Cu 1 -Cu 2 .
Hindwing: The subcostal ray on cell Sc+R 1 -Rs is expressed on the dorsum in full orange-red (expressed in pale whitish scales in H. e. ucayalensis). Also distinguishable from males by the five-segmented prothoracic tarsus (fused in male) and external genitalia.

Host plant ecology and immature morphology
In the wild near Tarapoto, confirmed host plant records for H. e. ucayalensis comprised clusters of 12-20 yellow ovoid eggs (N = 3), or groups of 1-4 gregarious larvae (N = 2) encountered on new leaves of Passiflora skiantha Huber (Passifloraceae: subgenus Astrophea) at Urahuasha (-6.466°, -76.335°) and San Roque de Cumbaza (-6.363°, -76.441°) (Fig. 5E, G). Both male and female H. e. ucayalensis were also often caught investigating P. skiantha plants in these and other nearby localities. When placed in an insectary with 22 local species of Passiflora, wild caught females (N = 6) laid 78 eggs on P. skiantha, in clusters of 12-33 eggs (N = 4), usually on new leaves and once on the expanding young shoot (Fig. 5A). One other female laid a single egg on Dilkea retusa Mast. (Passifloraceae). This latter female did also show considerable interest in P. skiantha prior to ovipositing on D. retusa, but the P. skiantha plant had no new growth at the time. Final instar larvae are characterized by a black head, legs and prolegs, spines and anal shield (Fig. 5B). Aside from the spiracles and a black band comprising a pair of elongated black spots running laterally on the dorsal side of the prothorax, only faint black spotting is observed on the thorax and abdomen, which are yellow. However, the larvae are notable for having black, annular stripes that start around the midpoints of each abdominal segment and run laterally and dorsally, approximately through the spiracles and the base of the spines. In between these black stripes, there are also fainter bands of darker coloration running between the abdominal segments. The pupae are typical for Heliconius in the H. erato clade, with long head horns (Fig. 5D). The base coloration is predominantly brown but with some paler bands/patches, and with distinct narrow white bands running horizontally and diagonally in the abdominal segments. There are three pairs of silver spots on the dorsal side of first abdominal segments, and an additional pair on the head. The horns are more darkly coloured and the spines are black. The horns are similar in length to those of H. erato and H. charithonia (Linnaeus, 1767), but are more elongate and taper to a point. Spines on the abdominal segments are somewhat longer than in H. erato and H. charithonia, and similar in length to those of H. sara.
Around Tarapoto we noted an association between presence of D. retusa and H. d. joroni ssp. nov. On several occasions H. d. joroni ssp. nov. females were caught in the vicinity of D. retusa plants at Biodiversidad (-6.460556°, -76.289928°), San Roque de Cumbaza, Pucayaquillo (-6.5882°, -76.2224°) and at La Antena (-6.45716°, -76.29858°). On two occasions, pairs of eggs were found on plants at Biodiversidad and La Antena; however, in general, finding eggs and larvae proved difficult. This is probably because it is difficult to find D. retusa with new growth suitable for the immature stages of Heliconius, at least in those plants accessible to human observers. Nonetheless, on 28 March 2016 a single first instar larva and a yellow ovoid egg were found on a D. retusa plant above San Roque de Cumbaza (Fig. 5F). The larva was reared to final instar using D. retusa (it refused P. skiantha), but failed to pupate. Its identity was confirmed as H. d. joroni ssp. nov. using COI DNA barcoding. This final instar larva was broadly similar to H. e. ucayalensis morphologically (Fig. 5C). However, the black annular stripes running between the spines were absent, and instead the larva was characterized by regular black spotting between the spines. The base colour also appeared a more greenish yellow than in H. e. ucayalensis; however, on the basis of a single individual it is unclear whether this is a reliable diagnostic character.
While we only provide data on larval morphology and host plant use from northern Peru, previously published data suggest that the specific differences we found in sympatry are widely applicable across the ranges of H. demeter and H. eratosignis. Heliconius demeter terrasanta Brown & Benson, 1975 has solitary, spotted final instar larvae and uses Dilkea sp. in the Brazilian state of Pará. Heliconius eratosignis eratosignis has been recorded using Passiflora ca. citrifolia Salisb. (subgenus Astrophea) in Rondônia, and has gregarious, striped final instar larvae (Brown & Benson, 1975).

Molecular data
The models of sequence evolution selected for each gene, along with associated parameter values and Bayesian Inference Criterion score, are shown in Supporting Information, Table S10. Analysis of mtDNA sequences (COI + COII) revealed a deep divergence between two haplogroups corresponding to H. d. joroni ssp. nov. and H. d. bouqueti + H. e. ucayalensis (Fig. 6). The net proportional distance between these haplogroups is 5.2%, and reciprocal monophyly was well supported (bootstrap percentages of 97% and 99%, respectively). Within the H. demeter cluster, H. d. bouqueti and H. d. joroni ssp. nov. also formed two well-supported, reciprocally monophyletic groups (bootstraps of 88% and 97%, respectively).
In addition, we were able to obtain ~760bp of COI sequence for 13 of the 15 previously recognized subspecies in the clade formed by H. demeter + H. eratosignis (Fig. 6). The resulting phylogeny indicated two reciprocally monophyletic groups, comprising the northern (H. demeter) and southern (H. eratosignis) races. The southern clade was well supported (100% bootstrap), and comprised H. e. ucayalensis, along with H. e. eratosignis, H. e. tambopata Lamas, 1985and H. e. ulysses Brown & Benson, 1975. The northern clade comprised H. d. demeter and H. d. bouqueti, along with H. d. angeli Neukirchen, 1997, H. d. karinae Neukirchen, 1990, H. d. neildi Neukirchen, 1997, H. d. terrasanta, H. d. titan and H. d. turneri Brown & Benson, 1975. The bootstrap support for this northern clade was only moderate (64%), however this is due to the uncertain placement of H. d. titan, which appears as sister to a well-supported (100%) monophyletic clade containing the other races of H. demeter.
Of the five nuclear loci examined, only Ef1α showed H. d. demeter + H. d. bouqueti and H. e. ucayalensis to form reciprocally monophyletic groups ( Supporting  Information, Fig. S6). Bootstrap support for these two groupings was only moderate (65% and 62%, respectively), and the two exhibited only two fixed nucleotide Between-run consistency was high in the STRUCTURE analysis of AFLP genotypes: replicate runs at each K-value yielded virtually identical likelihoods. The optimal number of genotypic clusters was three, corresponding cleanly to each of the three taxa (Fig. 7)

Geographic distribution
Subspecies of H. demeter and H. eratosignis are mapped in Fig. 1, with photos of a type specimen of each race. Races of H. demeter occupy the Guianas and much of the Amazon basin. H. eratosignis races occur in the west and south of the Amazon basin. In Tarapoto, the two species fly together at a number of sites in the Cordillera Escalera. Only H. eratosignis has been recorded from the adjacent Amazonian lowlands, despite considerable sampling in the area. Museum data and observations by Keith Brown (1979) suggest that the two overlap (at least broadly) in the extreme south of Pará and northern Mato Grosso, in Brazil. There may well also be a contact zone on the Juruá River, between Porto Walter and Eirunepé, as both H. demeter demeter and H. eratosignis tambopata are known to occur there. However, the exact position of contact in this very large area is unclear. In data published by Brown (1979) two additional contact zones are indicated, at Pucallpa, Peru and near Cobija on the Brazilian/Bolivian border. We were unable to locate the relevant specimens in museum collections; however, we consider these points unreliable and excluded them from the distribution map in Fig. 1. The first is probably a generalized locality, with the specimens potentially coming from a large area of northern Peru. The second is likely explained through the cooccurrence of both H. eratosignis ulysses and H. eratosignis tambopata, as the latter was not described at the time (Lamas, 1985).

DISCUSSION
Gene genealogies can be used in concert with morphological differences to diagnose species within single populations, because reciprocal monophyly within a freely interbreeding population becomes highly improbable when multiple individuals are sequenced. Similarly, the existence of clusters of multilocus genotypes within a sympatric population comprises strong evidence for distinct species, because linkage disequilibria between alleles at unlinked loci are highly unlikely to arise without barriers to recombination. We have shown that in northern Peru, H. d. joroni ssp. nov. and H. e. ucayalensis sampled from a small geographic area comprise two monophyletic groups for the mtDNA markers COI + COII, and form distinct genotypic clusters using AFLP data. Furthermore, the 5.2% net mtDNA divergence between H. d. joroni ssp. nov./H. d. bouqueti and H. e. ucayalensis is equivalent to interspecific genetic distances between other sarasapho group species, and is greater than distances between many other sister pairs of Heliconius species, such as those within the cydno-melpomene species group (Beltrán et al., 2002;Giraldo et al., 2008). Thus, together with the observed differences in larval and adult morphology, wing shape, behaviour and host plant use, our data strongly imply the existence of two species that are sympatric in at least one area. Additionally, the COI phylogeny of 13 of the 15 races of H. demeter and H. eratosignis resolved two reciprocally monophyletic groups, comprising H. d. joroni ssp. nov. and the northern Amazonian races, and H. e. ucayalensis and the southern Amazonian races. These groups are consistent with morphological criteria (e.g. the forewing costal streak) and are also evident in the morphometric analysis of wing shape ( Fig. 4; Supporting Information, Fig. S5). Both clades were well supported, excepting the uncertain position of H. d. titan, whose assignment to H. demeter rather than to H. eratosignis was only marginally favoured by molecular and morphological data. Heliconius demeter titan is also notable for discordant morphological characters, and for its long mtDNA branch lengths and reciprocal monophyly with H. demeter. Because H. d. titan appears broadly sympatric with other H. demeter races, it may even represent a further cryptic species within this clade.
In contrast to the mtDNA, only one of the five nuclear markers sequenced (Ef1α) showed reciprocal monophyly between H. d. bouqueti/H. d. joroni ssp. nov. and H. e. ucayalensis. However, two other nuclear genes (Tpi and Mpi) did show monophyletic groups corresponding to subspecies or species. Gene genealogies that fail to resolve relationships between closely related species are not unusual in Heliconius and may reflect either the retention of ancestral polymorphisms, introgression following speciation, or simply uninformative genetic data (Maddison, 1997;Beltrán et al., 2002;Bull et al., 2006). Because effective population sizes are lower for the maternally inherited COI + COII and sex-linked Tpi than for the autosomal loci, they are expected to coalesce more recently (Palumbi et al., 2001), and so finding monophyly at these loci remains consistent with the hypothesis of ancestral polymorphisms. In addition, if introgression was producing the observed patterns, we might expect polyphyly between the sympatric taxa H. e. ucayalensis and H. d. joroni ssp. nov., but with H. d. bouqueti phylogenetically distinct, due to its geographic isolation. However, females are the heterogametic sex in butterflies, and, in accordance with Haldane's rule, female sterility is an early manifestation of intrinsic postzygotic reproductive isolation (Jiggins et al., 2001a;Naisbit et al., 2002). Introgression should, therefore, be more inhibited at COI + COII and Tpi (Sperling, 1994), thus their monophyly could still be consistent with autosomal introgression between the species. As such, we cannot rule out introgression as a possible cause of incongruence between nuclear genealogies and species boundaries, especially given the abundant evidence for gene flow between closely related Heliconius (Dasmahapatra et al., 2007;Mallet et al., 2007;Heliconius Genome Consortium, 2012;Pardo-Díaz et al., 2012;Martin et al., 2013), and the known importance of colour pattern as a prezygotic reproductive isolating barrier in Heliconius (Merrill et al., 2011(Merrill et al., , 2012. Previous studies of cryptic Heliconius have suggested that hybridization between closely related co-mimics may be higher than between nonmimics, although quantitative comparisons are difficult (Giraldo et al., 2008;Mérot et al., 2013Mérot et al., , 2017. It would be interesting to investigate whether the other similarly divergent co-mimetic sister pair H. leucadia Bates, 1862 and H. sara exhibit similar phylogenetic discordance. It is striking that three recently described cryptic species pairs of Heliconius are distinguishable using a minor colour pattern difference in the costa-subcosta space on the forewing underside (Giraldo et al., 2008;Mérot et al., 2013;and the present study). Many other co-mimetic Heliconius are distinguishable using seemingly inconsequential red dots and streaks at the base of the ventral hindwing (Emsley, 1965;Holzinger & Holzinger, 1994). While this variation might be attributable to relaxed selection from predators on the underside of the hindwing, their repeated utility for distinguishing the species leads one to speculate that they are important for the butterflies themselves in terms of mate recognition. Indeed, these ventral areas are perhaps the most visible part of the wing to both sexes during courtship. This hypothesis could conceivably be tested using colour pattern manipulations and assortative mating experiments.
In other recently described cryptic Heliconius, phenotypic similarity is most parsimoniously explained by convergence through introgression of colour pattern alleles (Mallet, 2009;Heliconius Genome Consortium, 2012;Pardo-Díaz et al., 2012). In the case of H. demeter and H. eratosignis, the available data suggest that speciation occurred from start to finish without a significant mimicry shift. The present geographic distributions of the species are suggestive of vicariance between the north and south Amazon basin. This seems consistent with the species mimetic similarity, because allopatric speciation does not require ecological divergence (Coyne & Orr, 2004). It might also explain the poly-and paraphyly at nuclear loci, because monophyly would be slow to develop in the large vicariant populations (Maddison, 1997). Nonetheless, H. demeter and H. eratosignis do differ in other ecologically relevant traits that may have played a part in their speciation. Sexual dimorphism in colour pattern is very unusual in Heliconius, and finding that closely related species differ markedly in mating signals is often considered indicative of speciation via sexual selection (Panhuis et al., 2001). The 'greenish' scales (in reality, interspersed black and yellow scales) exhibited by males produce a seemingly non-mimetic phenotype that could be the product of sexual selection, but seem unlikely to be involved in speciation because they are present in both species. In contrast, fused rays are exhibited only by H. demeter. In some regions, such as near the Andes, this leads to males being somewhat poorer mimics of other Heliconius species than are females and could, therefore, be interpreted as the product of female choice for a male trait. However, in other regions, such as in French Guiana, the dimorphism seems to be a mixed strategy, with males mimicking species such as Heliconius egeria (Cramer, 1775) and females mimicking species such as H. erato. A mimetic explanation for the fused rays of H. demeter may, therefore, be more likely than sexual selection, and furthermore fits the hypothesis of vicariance, followed by more recent contact in the Amazon headwaters.
Heliconius demeter and H. eratosignis are also unusual in their apparent host plant specificity, because most Heliconius sister species use overlapping suites of Passiflora spp. (Rosser et al., 2015). Host plant shifts are frequently associated with speciation in phytophagous insects (Bush, 1969;Drès & Mallet, 2002), and there is some evidence for their importance in Heliconius (Jorge et al., 2011;Merrill et al., 2013;Rosser et al., 2015). This could be either because the butterflies tend to mate in the vicinity of their host plants (Bush, 1969), or due to disruptive selection for larval performance on alternative hosts (Funk, 1998). Heliconius demeter and H. eratosignis belong to a clade of Heliconius known to exhibit 'pupal mating', in which mating sometimes occurs on the host plant before the females have fully emerged from their pupae (Deinert et al., 1994), thus the former model seems possible. It also seems plausible that the evolutionary and phenotypic divergence between P. skiantha and D. retusa could produce disruptive selection on larval performance. For example, P. skiantha contains cyanogenic glycosides (secondary defence compounds) not found in D. retusa (Érika de Castro & Neil Rosser, unpublished). Furthermore, H. demeter and H. eratosignis are the only sister species pair within Heliconius known to comprise a species with gregarious larvae and one with solitary larvae (Beltrán et al., 2007;Kozak et al., 2015). Their larvae may also be involved in mimicry with other Heliconius species (Brown & Benson, 1975): Heliconius eratosignis larvae are nearly identical to the gregarious larvae of H. doris (Linnaeus, 1771) and H. xanthocles Bates, 1862(Brown & Benson, 1975Mallet & Jackson, 1980), whereas, H. demeter larvae are more similar to those of H. ricini (Linnaeus, 1758). Whatever the drivers of divergence in H. demeter and H. eratosignis, their limited geographic overlap, comimicry, sexual dimorphism, and marked differences in host plant use and oviposition behaviour, highlight them as an interesting counter-example to other Heliconius sister species. In particular, H. demeter and H. eratosignis exhibit striking parallels to cryptic species in the Afrotropical butterfly genus Cymothoe (Nymphalidae). Strong host plant and ecological differences have evolved between C. egesta (Cramer, 1775) and C. confusa Aurivillius, 1887, formerly considered subspecies of a single widely distributed species. These differences are apparently insufficient to allow sympatry, bar a narrow region of overlap between their otherwise allopatric ranges (McBride et al., 2009).
Our use of integrative taxonomy (Dayrat, 2005;Pante et al., 2015) to diagnose a cryptic species of Heliconius joins a series of similar, recent discoveries in other butterflies (Willmott et al., 2001;Hebert et al., 2004;McBride et al., 2009;Dincă et al., 2011;Hill et al., 2012;Barbosa et al., 2015). Frequently, cryptic taxa are initially flagged by molecular markers, after which subtle differences in morphology or behaviour are recognized as species-specific (Janzen et al., 2009). Thus, despite its limitations (Elias et al., 2007;Silva-Brandão et al., 2009), DNA barcoding still holds great potential to screen putative cryptic species for further study. While the net contribution of cryptic species to biodiversity remains to be established (Stork, 2018), the continual discovery of hidden species in a group as intensively studied as butterflies suggests that predictions of global species richness based on current knowledge may be gross underestimates (Adis, 1990;Bickford et al., 2007).

DATA ACCESIBILITY
DNA sequences have been submitted to GenBank; accession numbers are given in the Supplementary Information.

SUPPORTING INFORMATION
Additional Supporting Information may be found in the online version of this article at the publisher's web-site: Figure S1. Landmarks used for the analysis of wing shape and size with geometric morphometry (displayed on H. eratosignis ucayalensis male). Figure S2. Morphology of genitalia. Arrow indicates possible fixed difference. Voucher IDs in brackets. Figure S3. Forewing (left) and hindwing (right) centroid size for H. demeter joroni ssp. nov. and H. eratosignis ucayalensis. Groups labelled with the same letter do not show significant centroid size difference. Figure S4. Principal component analysis displaying wing shape variation between H. demeter joroni ssp. nov. (red) and H. eratosignis ucayalensis (blue). Males are represented with open circles and females by filled circles. Ellipses represent a graphical summary of the distribution. Top: forewing, bottom: hindwing. Shape variation is illustrated next to each axis, where broken shapes represent minimum negative values of the axis, and full lines represent maximum values. Figure S5. Principal component (PC) analysis of wing shape variation between H. eratosignis and H. demeter with the inclusion of additional subspecies from museum collections (types, syntypes, etc.). As in Fig. 4 H. d. angeli; b, H. d. bouqueti; d, H. d. demeter; k, H. d. karinae; n, H. d. neildi; t, H. d. titan; z, H. d. zikani; and with a blue letters for H. eratosignis: e, H. e. eratosignis; u, ucayalensis; y, H. e. ulysses. Ellipses represent a graphical summary of the distribution. Shape variation captured PC1 and PC2 are illustrated next to each axis, with dotted lines represent minimum values of the axis, and full lines representing the maximum values. PC1 captures shape differences between males and females across both species. PC2 captures variation between species as well as between H. demeter subspecies. Figure S6. Maximum likelihood trees of H. demeter and H. eratosignis specimens (red and blue respectively) from Tarapoto (Peru) and H. demeter from French Guiana (orange) based on sequences of five nuclear loci Mpi, Tpi, Tektin, Rpl5 and Ef1α. Bootstrap values greater than 50% are shown. Otherwise identical voucher numbers terminating in A or B refer to alleles from heterozygous individuals. Table S1. Revised synonymy of taxa formerly considered part of H. demeter. Letters a-m are valid subspecies names according to our revision and that of Lamas (2004).  .  Table S4. List of Passiflora species used in captive host plant oviposition tests. Table S5. Details of samples used for molecular work. Table S6. PCR conditions for the amplicons used in this study. All reactions were carried out in 10-µL volumes using 10×buffer (Sigma), 0.5 µM each of forward and reverse primers and 0.25U Taq polymerase (Sigma). PCR cycling conditions were: initial denaturing at 94 ºC for 2 min; followed by 35 cycles of 94 ºC for 45 s, annealing temperature (T a ) for 45 s, 72 ºC for 60 s; final extension at 72 ºC for 5 min. Table S7. GenBank Accession numbers for sequences used in nuclear phylogenies. Table S8. GenBank accession numbers for sequences used in mtDNA phylogenies. Table S9. Subspecies classified by presence or absence of strong costal yellow streak on the ventral forewing. Descriptions based on examination of holotype, syntype, allotype detailed in Table S2; fw = forewing; hw = hindwing. Table S10. Models of sequence evolution and estimated parameters, selected using Bayesian Inference Criterion.