Quantified reproductive isolation in Heliconius butterflies: Implications for introgression and hybrid speciation

Abstract Heliconius butterflies have become a model for the study of speciation with gene flow. For adaptive introgression to take place, there must be incomplete barriers to gene exchange that allow interspecific hybridization and multiple generations of backcrossing. The recent publication of estimates of individual components of reproductive isolation between several species of butterflies in the Heliconius melpomene–H. cydno clade allowed us to calculate total reproductive isolation estimates for these species. According to these estimates, the butterflies are not as promiscuous as has been implied. Differences between species are maintained by intrinsic mechanisms, while reproductive isolation of geographical races within species is mainly due to allopatry. We discuss the implications of this strong isolation for basic aspects of the hybrid speciation with introgression hypothesis.


| INTRODUCTION
Since the demise of Ernst Mayr (1904Mayr ( -2005, emphasis on allopatry in the study of speciation has been eclipsed by ideas about ecological factors that might drive divergence of taxa in the absence of extrinsic barriers to reproduction (Jiggins, 2008;Rundle & Nosil, 2005;Schluter, 2009). The geographical component of speciation, once viewed as a key factor permitting the initial stages of divergence, which would otherwise be swamped by gene flow, is now often viewed as unnecessary, or at least passé. Various animals, including sticklebacks, timemas, and Darwin's finches, have become textbook icons exemplifying this shift of theoretical and empirical focus (Nosil, 2012).
Heliconius butterflies have likewise become a model system for studying patterns of speciation in the putative presence of ongoing gene flow Kronforst et al., 2006Kronforst et al., , 2013Nadeau et al., 2013). Most Heliconius species are Müllerian mimics of one another, as well as other unpalatable taxa, displaying shared, aposematic wing patterns that are maintained by positive numerically dependent selection. Although selection would seem to favor a single, widespread signal to potential predators, paradoxically, multiple mimetic patterns may exist among different species occurring at a given locale, and further, many of these species also exhibit dramatic geographical variation in these wing patterns, so that the particular pattern shared between butterflies that mimic one another varies from place to place across their Neotropical distributions (Turner, 1975). While positive numerically dependent selection maintains this variability, it cannot explain its origin.
Traditionally, mimetic resemblance has been explained by convergent evolution, with similar wing patterns arising independently in separate lineages. It is clear, based on the diversity and phylogenetic distribution of particular aposematic patterns, that convergence remains the most plausible explanation for many instances of mimicry both within Heliconius and between the genus and its comimics ( Figure 1). Thus, although convergence is a less parsimonious explanation for similarity than common ancestry, when mimetic wing patterns have evolved independently on multiple occasions within a taxon (not to mention in more remotely related groups- Figure 1), it is less onerous to explain shared patterns by that common mechanism than it is by invoking further, novel alternatives. One ad hoc hypothesis is more parsimonious than two ad hoc hypotheses.
Despite this, recent research has offered a new explanation for mimetic resemblance, suggesting that some Heliconius populations may have been able to shift from one mimicry ring to another by virtue of hybridization leading to adaptive introgression of mimetic wing pattern alleles and that one or more taxa may have arisen as a result of homoploid hybrid speciation Enciso-Romero et al., 2017;Mávarez et al., 2006;Pardo-Diaz et al., 2012).
These hypotheses are counterintuitive, given the aforementioned F I G U R E 1 Multiple origins of the "dennis-ray" mimetic pattern in Heliconius and other Lepidoptera. Cladogram of Heliconius species based on Brower and Garzón-Orduña (2017). Exemplar Heliconius exhibiting the dennis-ray pattern are illustrated (top to bottom, H. erato, H. demeter, H. aoede, H. doris, H. burneyi, H. melpomene, H. timareta, H. elevatus), and the origins of those features are parsimoniously optimized, indicating at least eight separate origins (red branches on tree and taxon labels). Note that H. erato, H. melpomene, and H. timareta also include geographical races that do not exhibit the dennis-ray pattern. Exemplars of Müllerian or Batesian mimetic phenotypes of more distantly related butterflies and moths representing at least six further independent origins of the dennis and/or ray pattern are inset. Images are open access (Wallbank et al., 2016) or courtesy of Keith Willmott, Florida State Museum of Natural History selective regime that is thought to establish and maintain mimetic patterns in Heliconius communities. If positive numerically dependent selection favors abundant aposematic patterns, then if they could do so, why would not all potentially hybridizing species (not to mention actually hybridizing geographical races of the same species) converge upon the same mimetic pattern? One reason that wing pattern diversity is maintained could be that interspecific hybridization in Heliconius may not be as common and straightforward as has been advanced in the literature (Mallet, Beltrán, Neukirchen, & Linares, 2007).
According to the most recent checklist (Lamas & Jiggins, 2017), the melpomene-cydno clade (Figure 1) comprises five species, three of which exhibit extensive diversification into phenotypically differentiated geographical races (Arias et al., 2014;Brower, 1996;Brown, 1979 the Andes, where it is now recognized to occur in at least seven geographically differentiated races that for the most part exhibit red and yellow pattern elements that mimic to a greater or lesser degree the sympatric races of H. melpomene and H. erato (Brower, 2011(Brower, , 2013Mallet, 2009). Nested allopatrically among these is H. heurippa, which occurs on the eastern slopes of the Andes in central Colombia, and although it has red and yellow forewing bands, is considered to be nonmimetic, as its wing pattern is not like that of any other Heliconius species. Although a consensus is beginning to develop that H. heurippa and H. timareta forms are conspecific (Arias et al., 2014;Brower, 2011;Jiggins, 2017; in which case all H. timareta forms should be considered subspecies of H. heurippa due to nomenclatural priority of publication), the traditional species names will be employed here.
Basic prerequisites for exchange of alleles between species are the success of interspecific mating events that produce F1 hybrids and the subsequent mating success of those hybrids with one or the other of the parental forms, resulting in introgression of alleles from one parental species into the other (Rieseberg & Wendel, 1993). It has long been known that mate discrimination in Heliconius is facilitated by visual recognition cues based on wing patterns (Crane, 1955;Jiggins, Naisbit, Coe, & Mallet, 2001;Merrill et al., 2011). This serves not only to promote conspecific, and to deter heterospecific mating in parental forms, but also to preclude mating opportunities for hybrid offspring . Further, interspecific mating experiments have shown that closely related species often produce sterile female offspring (Nijhout, Wray, & Gilbert, 1990;Salazar et al., 2004) following Haldane's Rule. Beyond these reproductive constraints, hybrid offspring with novel combinations of wing pattern elements that disrupt their participation in established mimicry rings also suffer strong selection due to predation (Mallet & Barton, 1989). All of these phenomena would seem to strengthen selective barriers to interspecific gene exchange (Brower, 2011).

| METHODS
We used Mérot et al.'s (2017) estimates of individual isolation components (their Table 1) to calculate total reproductive isolation (TI) using the methods proposed by Sobel and Chen (2014). Briefly, Sobel and Chen suggested describing the relationship between the probability of gene flow and the probability of reproductive isolation as a linear equation that expresses the probability of reproductive isolation from 0 to 1, the former representing unrestricted gene flow (or disassortative mating) and the latter complete reproductive isolation (probability of gene flow under random mating is .5). Currently, data on the various components of reproductive isolation are available only for 12 matings involving various geographical races of the species described above; the comparisons examined here are based on these 12 mating pairs (Figure 2). The names and localities of these comparisons are given in Table 1. These comparisons range from interspecific (pairs 1-8), to between geographical races (pairs 9-11), and to sympatric polymorphic forms (pair 12).
Sobel & Chen's calculation categorizes barriers to gene flow into three types: (extrinsic) prezygotic barriers that affect co-occurrence, (intrinsic) prezygotic barriers not related to co-occurrence, and postzygotic barriers (Sobel & Chen, 2014). As applied to the Heliconius species tested here, these variables reflect the degree of spatial overlap of the parental species, habitat preference, interspecific mating success, and the viability, mating success, and fertility of hybrid offspring. The values for the spatial component and mating were assigned to Sobel & Chen's first and second categories, respectively, while all the other variables were considered postzygotic barriers (F1 adult, F1 fertility, F1 egg, F1 mating with parent # 1, F1 mating with parent #2). We note, however, that the spatial co-occurrence values of Mérot et al.
(2017) reflected a conflation of allopatry in the traditional (geographical) sense and ecological habitat preference, such that species that appear to be sympatric according to range maps might have a deficit of encounters with each other due to relatively small altitudinal or ecological differences, such as larval host plant preference. For example, they reported a spatial isolation value of 0.74 between H. melpomene rosina and H. cydno chioneus, even though samples of both were collected along a few kilometers of Pipeline Road in Panama. As Heliconius butterflies are vagile animals that may move several kilometers over the course of their lives (Mallet et al., 1990), we view habitat preference within a geographical area as an intrinsic rather than an extrinsic barrier, and have separated these components accordingly in Table 2. Mérot et al. (2017) provided values for the strength of spatial isolation for only four of 12 interspecific mating comparisons, and we complemented these for the other eight species pairs based on distributional data regarding sympatry, parapatry, or allopatry (Brown, 1979 Phenotypes are illustrated in Figure 2. a In some instances, stocks for a given comparison were founded from specimens collected at more than one site. See the cited references for details. b Data reported incorrectly by Mérot et al. (2017). component, which reflects the "actual" and "potential" aspects of reproductive isolation of the Biological Species Concept (Mayr, 1942). We emphasize that these values should be viewed as probabilities of reproductive isolation, not the inverse frequency or rate of interbreeding.

| RESULTS AND DISCUSSION
Estimates of total reproductive isolation in Table 2 (Kronforst, Young, & Gilbert, 2007;Kronforst et al., 2006). In contrast, H. heurippa is only isolated from H. cydno by allopatry, suggesting that these two forms are parts of a single biological species (Brower, 2011 (Kronforst et al., 2007).
Although much has been made of the importance of color pattern in Heliconius mate recognition Mávarez et al., 2006), the values for total isolation in Table 2 show that the two comimetic and visually almost identical, sympatric pairs of H. melpomene and H. timareta (pairs 7 and 8) are among the most isolated taxa in the dataset, and that a major contributing factor to this isolation is mate choice. This indicates that there must be additional cues, such as behavioral or pheromonal differences that allow species recognition/discrimination even when the butterflies look virtually the same (Mérot, Frérot, Leppik, & Joron, 2015). Nonvisual stimuli could serve both to enhance probabilities of intraspecific mating and to deter interspecific mating (Darragh et al., 2017;Friberg et al., 2008).
Experiments testing only the visual component of mate choice using paper wing models, or using freshly emerged virgin females that may not have a full behavioral repertoire, are therefore not only artificial, but do not present a complete assessment of potential components of mate choice and/or discrimination. Thus, the "total isolation" values shown in Table 2 are likely to be underestimates of the actual isolation, when differential chemistry and behavior contribute to mate choice in nature. A further observation supporting the role of nonvisual cues in mate recognition is the fact that numerous hybrid zones exist within both H. melpomene and H. cydno, in which phenotypically different geographical races of the same species freely interbreed despite the fact that they have different wing patterns.
Of course, experimental evidence of "complete" reproductive isolation does not preclude the possibility of rare interspecific mating events. Through evolutionary time, extremely unlikely events can play a significant role. However, even if hybridization occurred occasionally, putatively beneficial introgressed alleles would need to survive several additional rounds of backcrossing to become integrated into the opposite genome. F1 and subsequent backcross offspring could suffer comparable pre-and postzygotic losses of fitness due to mate discrimination, sterility, and predation as measured above. If that is true, the quantitative values in Table 2 thus represent an overestimate of the very small likelihood of introgression.
The hypothesis stated by Mérot et al. (2017), "reproductive isolation between pairs at a high level of divergence is strong enough to allow the secondary loss of certain barriers to gene flow, in this case via the introgression of wing pattern alleles, without compromising genome-wide differentiation," may provide a plausible argument for the maintenance of mimicry, but not for its origin. By definition, if a species gains a new wing pattern via introgression, then it must have started with a different wing pattern than the one it displays today.
The loss of the visual component of preexisting barriers to gene flow is a consequence of mimetic convergence, but necessitates that other barriers are strong enough to maintain species integrity, lest the two taxa undergo reverse speciation (Lackey & Boughman, 2016), and further implies that the species were even more strongly isolated from one another prior to the gene flow than they are now.
It has been suggested that H. heurippa and H. timareta forms have obtained their red pattern elements by introgression from H. melpomene (Giraldo, Salazar, Jiggins, Bermingham, & Linares, 2008;Mávarez et al., 2006;Pardo-Diaz et al., 2012). which H. timareta has acquired its red wing pattern elements, then that mechanism fails to explain why H. timareta timareta is polymorphic with nonmimetic forms, and for that matter, why H. heurippa is not a comimic of its sympatric H. melpomene race (Brower, 2011).
If interspecific introgression of wing pattern alleles is a real phenomenon in Heliconius, then the interplay of selection and gene flow must be fundamentally different in cases of intraspecific versus interspecific hybridization. In intraspecific hybrid zones (e.g., the H. melpomene amaryllis-H. melpomene aglaope zone studied by Mallet & Barton, 1989; in the Huallaga Valley of Peru, Pair 11 in Table 2), there is no apparent mate discrimination based on wing pattern, and interracial mating takes place quite freely (Figure 3a). The two geographical races are maintained as distinct by very strong positive numerically dependent selection acting on those alleles, while alleles not related to wing pattern apparently mix readily (Turner, Johnson, & Eanes, 1979), such that there are "islands of divergence" in a sea of genetic homogeneity . In contrast, in interspecific crosses (Figure 3b), gene flow is minimal, due to the rarity of hybridization events, and wing pattern alleles become islands of similarity in a genomic sea that reflects underlying phylogenetic relationships (Enciso-Romero et al., 2017). The apparently contradic- Traditionally, the origin of Heliconius intraspecific wing pattern diversity has been explained by vicariance, perhaps due to the fragmentation of the rainforest during cool, dry Pleistocene climate cycles (Brown, 1979;Brown, Sheppard, & Turner, 1974). Molecular clock estimates for diversification of geographical races are consistent with this time frame (Brower, 1994;Garzón-Orduña, Benetti-Longhini, & Brower, 2014). Under this scenario, phenotypes of smaller, isolated populations could evolve novel wing patterns due to genetic drift (as in Phase 1 of the shifting balance; Wright, 1977), or be selected to converge upon wing patterns of other locally abundant unpalatable butterflies, such as members of the genera Altinote, Melinaea, or Elzunia (Turner & Mallet, 1996). Given the fact that all the current mimicry rings appear to have arisen relatively recently, it is possible that many other mimetic patterns may have existed in the past and gone extinct (cf. Linares, 1997). The strength of reproductive isolation found between sympatric species versus that between races separated only by geography is consistent with the hypothesis that the origin of novel wing patterns now maintained by intrinsic barriers was facilitated by allopatry.
Finally, we note that although the data on reproductive isolation reported by Mérot et al. (2017) and discussed here hardly characterize current species boundaries in the melpomene-cydno clade as a "continuum" (cf. Mallet et al., 2007), several objections might be posed.
First, the mating experiments might measure unnatural behavior, or be insensitive to rare events that could allow gene flow to take place in spite of empirical "total isolation." Or maybe the Sobel-Chen equation overestimates isolation. Third, perhaps contemporary strong reproductive isolation does not reflect the degree of isolation that may have existed over the past few hundred thousand years when purported T A B L E 2 Components of reproductive isolation and total isolation as calculated by the formula of Sobel and Chen (2014)

CONFLICT OF INTEREST
None declared.

AUTHOR CONTRIBUTIONS
IJGO and AVZB conceived the study, wrote the manuscript, and drafted the figures. IJGO ran the analyses.