\chapter{\label{cha:hybrids}Case studies and mathematical models of ecological
speciation.\protect \\
Hybrid speciation in butterflies in a jungle}
This chapter is a lightly revised version of a paper by the same name
to be published in the \emph{International Journal of Organic Evolution}
in 2009 by Edgar A. Duenez-Guzman, Jes\'us Mav\'arez, Michael D.
Vose, and Sergey Gavrilets:
Duenez-Guzman, E.A., Mavarez, J., Vose, M.D., and Gavrilets, S. Case
studies and mathematical models of ecological speciation. 4. Hybrid
speciation in butterflies in a jungle. \emph{Evolution}, 2009, in
press.
My use of \textquotedbl{}we\textquotedbl{} in this chapter refers
to my co-authors and myself. My primary contributions to this paper
include (1) discussion on the setup of the model, (2) implementation
of the numerical simulation and secondary analysis tools, (3) data
acquisition and analysis, (4) most of the data interpretation and
discussion, and (5) most of the writing.
\section{Introduction}
Our understanding of the processes leading to the origin of new species
has been steadily increasing both from empirical and theoretical perspectives
(e.g. \citealt{coy04,gav04,die04,mal07,nos07,see06,see07,bol07,gav09a}).
One of the lessons of recent work is a renewed appreciation of the
role of ecology in speciation which has lead to a new focus on \emph{ecological
speciation} (e.g. \citet{may47,schl00,run05}), that is, speciation
driven by ecologically-based divergent selection. Selection is divergent
when opposite phenotypes are favored within different populations
or the same population. Selection is ecological when it arises as
a consequence of the interaction of individuals with their abiotic
and biotic environment. Another novel development is a resurrection
of arguments about the role of hybridization in speciation and adaptive
radiation \citep{mal07,gom06,arn97,bul94,see04}. In particular, the
role of homoploid \emph{hybrid speciation}, i.e. hybrid speciation
without change in ploidy level (also referred to as \emph{recombinational
speciation}, \citealt{gra85}), is an issue receiving renewed interest
and new empirical and theoretical support \citep{mav06,nic02,gro05,gom06,sal04,see04,mav08}.
An important question is what is really meant by hybrid speciation.
We will call the result of introgressive hybridization a hybrid species
if resulting hybrid traits directly and significantly contribute to
the survival and reproductive isolation of the species \citep{mal07,mal09}.
Note that we require neither that genomes of the parental species
are represented in the hybrid species at equal frequencies nor that
the hybrid species stably coexist in sympatry with the parental species.
We find this relaxed view of hybrid speciation more useful than alternatives
as it fits better the patterns observed in nature. Indeed most suggested
cases of homoploid hybrid speciation among animals appear to have
involved a certain degree of backcrossing to one parental species
\citep{mav08}. Moreover these authors also noted that about half
of the homoploid hybrid taxa considered are allopatric with at least
one parental species.
Although the support for the importance of ecological speciation and
homoploid hybrid speciation in nature is growing, many questions remain.
These concern the conditions for speciation, its time-scales, driving
forces, the importance of ecological and genetic details, the role
of geography and so on. Answering these questions requires more detailed
data on already studied cases, new empirical studies, and extensions
of the quantitative theory of speciation.
One recent trend in theoretical research is to use complex numerical
simulations tailored to particular case studies to supplement and
provide additional insights to those that have emerged from relatively
simple analytical models \citep{gav04}. For example, recently \citet{gav07a,gav07b,sad08}
used detailed simulation models incorporating relevant ecological,
behavioral, spatial and genetic data to examine putative cases of
ecological speciation of cichlids in a crater lake \citep{bar06a},
of palms on an oceanic island \citep{sav06}, and of snails on sea
shores \citep{hol05,hol06}. By simulating these real systems, they
were able to address certain questions about ecological speciation
in general (e.g. whether sympatric speciation is achieved easily as
it is claimed occasionally) and in particular case studies (e.g. whether
an observed pattern is a result of \emph{in situ} speciation or double
invasion, whether speciation was truly sympatric or parapatric). They
also looked at specific questions such as: what does mathematical
theory tell us about the plausibility, speed, and patterns of (sympatric)
speciation in the case studies? What are the important parameters
and processes controlling the dynamics of speciation? How common are
the phenomena observed in these case studies? They also identified
some important parameters and features that need to be studied empirically
to provide information that can be used to improve the biological
realism and power of mathematical models of ecological speciation
and to make the interpretation of empirical findings less speculative.
In contrast to a significant effort to understand ecological speciation,
recombinational speciation has so far received only very limited attention
from theoreticians. Two previous theoretical papers by \citet{mcc95}
and \citet{bue00} have utilized mathematical models to help understand
the dynamics of hybrid speciation in plants. \citet{mcc95} used spatially
explicit individually-based simulations to study the consequences
of hybridization between two parental forms that differed in two chromosomal
rearrangements. These authors assumed that one homozygous hybrid genotype
({}``hybrid species'') had higher fertility whereas all heterozygous
hybrid genotypes had smaller fertility than both parental genotypes.
\citet{mcc95} analyzed how the waiting time until the advantageous
recombinant type is fixed in the population depended on parameters
(such as the strength of selection, the rate of selfing, and the size
of the area where hybridization occurs). \citet{mcc95} noted that
even when the relative fertility of $F_{1}$ hybrids was very low
(e.g. less than 2\%), the new species was established in just a few
hundred generations. Decreasing the hybrid species advantage markedly
increased the amount of selfing required for rapid speciation, but
when the advantage was sufficiently large, speciation was observed
even for obligate outcrossers. Increasing the size of the area where
hybridization occurs significantly reduced the time to speciation.
In the simulations of \citet{mcc95}, the environment was spatially
homogeneous and the new recombinant species, once emerged, always
replaced the parental forms. \citet{bue00} considered a more complex
situation with a spatially heterogeneous environment and an additional
{}``ecological'' trait controlled by two unlinked additive diallelic
loci. Specifically, they assumed a gradient in viability selection
so that small values of the ecological trait (which were characteristic
of one parental form) were advantageous in one habitat whereas large
values of the trait (which were the characteristic of another parental
form) were advantageous in the other habitat. The intermediate values
of the ecological trait resulted in higher viability in a relatively
narrow intermediate zone positioned between the parental habitats.
\citet{bue00} showed that under certain conditions the hybrid species
can appear and be stably maintained in the intermediate zone simultaneously
with the parental forms being preserved in their respective habitats.
The time scale for speciation was on the order of several hundred
to several thousand generations. An analysis of neutral markers showed
a significant gene flow and loss of differentiation in neutral loci
between the species.
Here, we build on the previous theoretical studies of hybrid speciation
\citep{mcc95,bue00} and ecological speciation \citep{gav07a,gav07b,sad08}
to analyze a tentative case of ecological hybrid speciation in \emph{Heliconius}
butterflies \citep{mav06,sal05}. This case has received a lot of
attention after the publication of a paper by \citet{mav06} arguing
in favor of hybrid speciation of \emph{H. heurippa} from \emph{H.
melpomene} and \emph{H. cydno} in South America. Our model attempts
to account for empirical patterns and data on genetic incompatibility,
mating preferences and selection by predation (both based on coloration
patterns), habitat preference, and local adaptation for all three
\emph{Heliconius} species. Using this model we study the likelihood
of recombinational speciation and identify the effects of various
ecological and genetic parameters on the dynamics, patterns, and consequences
of hybrid ecological speciation.
\section{Empirical evidence}
In describing empirical evidence we will closely follow \citet{mav06}
and references therein. \emph{Heliconius cydno} and \emph{H. melpomene}
are two closely related species that overlap extensively in lower
Mesoamerica and the Northern range of the Andes (i.e. Venezuela, Colombia
and Ecuador). Speciation in these butterflies does not involve any
change in chromosome number \citep{bro92} but it is instead associated
with shifts in wing color patterns that generate assortative mating
as well as post zygotic isolation due to predation-mediated selection
\citep{jig04,mal98,mcm97}. Different geographic races of \emph{H.
cydno} have black wings with white and yellow marks, whereas races
of \emph{H. melpomene} have black wings with red, yellow and orange
marks (see Figure \ref{fig:wings}). %
\begin{figure}
\hfill{}\includegraphics[width=6in]{Figures/hybrid1}\hfill{}
\caption[~~~The wing patterns of the three \emph{Heliconius} butterfly species]{\label{fig:wings} The wing patterns of the three \emph{Heliconius}
butterfly species. Some races of \emph{H. cydno} have black wings
with white and yellow marks in their forewings and a brown pincer-shaped
mark on their hind-wings (not visible). Races of \emph{H. melpomene}
have black wings with red and orange marks in their fore-wings and
brown dots in their hind wings (not visible). \emph{H. heurippa},
shown in the middle, has both yellow and red bands in their fore-wings,
and brown dots in their hind-wings (not visible).}
\end{figure}
Both species exhibit strong positive assortative mating based on
their wing patterns \citep{jig01,mav06} and also differ in habitat
use (\emph{H. cydno} in forest understory, \emph{H. melpomene} in
forest gaps/open areas) and the degree of preference for host plants
in the genus \emph{Passiflora} (low preference in \emph{H. cydno},
high preference in \emph{H. melpomene}). In spite of these differences,
interspecific hybrids still occur in the wild throughout their sympatric
range, usually at low frequencies ($<10^{-3}$, \citealt{mal98})
but sometimes representing a significant fraction of the population
($\sim10\%$, \citealt{mav06}). The tentative hybrid species, \emph{H.
heurippa}, is ecologically most similar to \emph{H. cydno}, which
it replaces geographically in the eastern Andes of Colombia. The wing
pattern of \emph{Heliconius heurippa} has elements of those of the
geographically close races of \emph{H. cydno} and \emph{H. melpomene}.
Its hind-wing is indistinguishable from that of sympatric \emph{H.
m. melpomene}, whereas its forewing shows a mixture of elements of
both \emph{H. m. melpomene} and of parapatric \emph{H. cydno cordula}.
Genetic evidence for a hybrid origin of \emph{H. heurippa} comes from
analysis of polymorphism at the two nuclear genes Invected and \emph{Distal-less}.
These genes show no allele sharing between \emph{H. cydno} and \emph{H.
melpomene}, whereas the \emph{H. heurippa} genome appears as an admixture,
sharing allelic variation from both putative parental species. Moreover,
by performing inter-specific crosses between \emph{H. cydno} and \emph{H.
melpomene} in the lab, \citet{mav06} were able to recover a color
pattern phenotype almost indistinguishable from that of wild \emph{H.
heurippa} in three generations, which provides some insight into the
possible steps of introgressive hybridization that could have given
rise to this species. Furthermore, the laboratory reconstructed pattern
breeds true when crossed among themselves and with wild \emph{H. heurippa}
individuals.
Most interspecific crosses in \emph{Heliconius} follow Haldane's rule
for sterility, i.e. $F_{1}$ females, which are the heterogametic
sex in Lepidoptera, are sterile. On one hand, crosses between \emph{H.
cydno} and \emph{H. melpomene} produce $F_{1}$ female sterility in
both directions of the cross. On the other hand, $F_{1}$ female hybrids
between \emph{H. heurippa} and \emph{H. melpomene} show asymmetrical
sterility. Finally, The female offspring of a cross between a female
\emph{H. heurippa} and male \emph{H. melpomene} are completely sterile,
either failing to lay eggs or laying eggs that never hatch. In contrast,
the reciprocal cross produces female offspring that lay fertile eggs.
Finally, crosses between \emph{H. cydno} and \emph{H. heurippa} produce
fertile $F_{1}$ females in both directions. These results are also
compatible with the hypothesis of hybrid origin for \emph{H. heurippa}
\citep{sal05}. In addition, back-crosses using $F_{1}$ hybrid males
provide evidence for a large Z(X)-chromosome effect on sterility and
for recessive autosomal sterility factors as predicted by the dominance
theory \citep{jig01,nai02,sal05}.
\citet{mav06} also studied genetic isolation on several sympatric
and allopatric populations of the three species using Bayesian assignment
analysis performed with 12 microsatellite loci. The three species
clearly represent different genetic entities, even in sympatry, which
suggests the operation of some form of reproductive isolation. To
further explore this idea, \citet{mav06} tested the degree to which
\emph{H. heurippa} is isolated from \emph{H. melpomene} and \emph{H.
cydno} by assortative mating. No-choice mating experiments were conducted
and showed that both males and females courted their co-specifics
up to ten times more often than individuals from other species. In
mating experiments with choice, there was similarly strong assortative
mating, although occasional matings between \emph{H. heurippa} and
\emph{H. cydno} were observed.
\section{Model description}
The following describes the major components of our model.
\textbf{Space and environment.}\quad{}Space is subdivided into a
$S_{x}\times S_{y}$ rectangular array of {}``patches'' each representing
roughly one square kilometer of forest. In Central America, \emph{H.
melpomene} are found in open habitat, \emph{H. cydno} are found in
closed-canopy forest \citep{smi78,est02}, and \emph{H. heurippa}
uses the same habitat as \emph{H. cydno} (Mavarez, pers.observation).
In our model, the environment changes in the horizontal direction.
We assume that both habitats (denoted as $0$ and $1$) are present
in the $S$ central columns of patches, whereas the $(S_{x}-S)/2$
leftmost columns and $(S_{x}-S)/2$ rightmost columns of patches have
only habitat $0$ ({}``open'') and habitat $1$ ({}``closed-canopy''),
respectively (see Fig~\ref{space}). Each habitat has a number of
host plants that can be used for oviposition. Because the central
area has more host plants present, the population density there will
be higher as well (see below). To reduce boundary effects, we wrap
the rectangle into a tube so that the top and the bottom rows become
neighbors. Time is discrete and generations are non-overlapping.
%
\begin{figure}[t]
\hfill{}\includegraphics[width=2in]{Figures/systemBW}\hfill{}
\caption[~~~Schematic representation of the spatial structure of the system]{\emph{ }Schematic representation\emph{ }of the spatial structure
of the system. Each rectangle represents a patch. Open circles represent
hosts in the open habitat; filled circles represent hosts in the closed-canopy
habitat. Note the presence of both types of hosts in the central area. }
\label{space}
\end{figure}
\textbf{Individuals.}\quad{}Individuals are sexual and diploid and
have discrete sexes. Each individual has a number of various phenotypic
characters. To capture the wing color pattern differences between
\emph{H. melpomene}, \emph{H. cydno}, and \emph{H. heurippa} we assume,
following \citet{mav06}, two unlinked diallelic loci (with alleles
$A,a$ and $B,b$) controlling the presence and intensity of red and
yellow bands on the dorsal fore-wing. Let $i$ and $j$ be the numbers
of {}``red'' ($A$) and {}``yellow'' ($B$) alleles in the genotype
($i,j=0,1$, or $2$). Then the intensities of red and yellow colors
are characterized by color intensity matrices \[
R=\left(\begin{array}{c|ccc}
& bb & Bb & BB\\
\hline aa & 0 & 0 & 0\\
Aa & 1 & 1 & 1/2\\
AA & 1 & 1 & 1\end{array}\right),\quad Y=\left(\begin{array}{c|ccc}
& bb & Bb & BB\\
\hline aa & 0 & 1/2 & 1\\
Aa & 0 & 0 & 1\\
AA & 0 & 0 & 1\end{array}\right).\]
These matrices capture epistatic interactions between the loci. The
two ones in the middle row of the $R$-matrix account for strong expression
of red color in some $Aa$ heterozygotes. Specifically, red allele
$A$ behaves as dominant in $bb$ and $Bb$ individuals, but its effect
are additive in $BB$ individuals. The two zeroes in the middle column
of the $Y$-matrix account for the absence of yellow color in some
$Bb$ heterozygotes (notice the lack of yellow band in the right box
in Figure 2a in the Mavárez et al. paper). Specifically, yellow allele
$B$ behaves as recessive in $Aa$ and $AA$ individuals, but its
effect is additive in $aa$ individuals. Under this parametrization,
\emph{H. melpomene} is represented by genotype $(i=2,j=0)$ with $R_{20}=1,Y_{20}=0$
(bottom left corner of the color intensity matrices), \emph{H. cydno}
by genotype $(0,2)$ with $R_{02}=0,Y_{02}=1$ (top right corner),
and \emph{H. heurippa} by genotype $(2,2)$ with $R_{22}=1,Y_{22}=1$
(bottom right corner). We note that there are several known color
pattern genes that differentiate \emph{H. melpomene} and \emph{H.
cydno}; we focus on only two of them for simplicity.
In Lepidoptera, males are the homogametic sex (and have two $Z$ chromosomes)
while females are the heterogametic sex (and have one $Z$ and one
$W$ chromosome). To account for Haldane's rule patterns observed
among the three species (\citealt{sal05}; see Introduction), we assume
that hybrid female sterility is due to incompatibility between $Z$-linked
and autosomal genes \citep{orr97}. Specifically, we postulate that
there are three types of sex chromosome $Z$ ($Z_{1},Z_{2},Z_{3}$)
and a diallelic autosomal locus with alleles $C_{1}$ and $C_{2}$
such that relative female fertilities are given by a matrix
\[
\begin{array}{c|ccc}
& C_{1}C_{1} & C_{1}C_{2} & C_{2}C_{2}\\
\hline Z_{1} & 1 & 0 & 0\\
Z_{2} & 1 & 1 & 1\\
Z_{3} & 0 & 0 & 1\end{array}.\]
Since the matrix concerns females exclusively, the $W$ sex chromosome
is omitted. \emph{H. melpomene} and \emph{H. cydno} females are represented
by the top left and bottom right corners, respectively. We assume
there is a small remaining polymorphism on the $Z$ chromosome to
account for incomplete divergence between \emph{H. cydno }and \emph{H.
melpomene}.\emph{ H. heurippa} females are represented by genotype
$Z_{2}C_{2}C_{2}$. This matrix simultaneously captures the bi-directional
Haldane's rule observed in hybrids between \emph{H. melpomene} and
\emph{H. cydno} (so that all $F_{1}$ hybrid females are sterile)
and the one-directional Haldane's rule observed in hybrids between
\emph{H. melpomene} and \emph{H. heurippa} (so that only the $F_{1}$
female offspring of male \emph{H. melpomene} $\times$ female \emph{H.
heurippa} crossings are sterile) \citep{mav06,sal05}. To keep the
model's complexity at a reasonable level, we neglect a possibility
that fertilities of fertile females differ and that there may be more
autosomal loci involved in reproductive isolation. All males are assumed
to be fully fertile. We note that although there seems to be some
sterility in males from certain hybrids crosses, it is definitely
not as strong as in females.
To account for mating preferences we assume the existence of two additive
quantitative characters $p_{r}$ and $p_{y}$ controlling males' preference
for the presence in females of the red and yellow fore-wing bands,
respectively. To allow for habitat preference and local adaptation,
we assume the existence of two additive quantitative characters: a
{}``habitat preference'' character $h$ controlling the probability
of choosing one or another habitat and an {}``ecological'' character
$x$ controlling the adaptation to host. The traits $p_{r},p_{y},h,$
and $x$ are scaled between $0$ and $1$ and are controlled by different
unlinked diallelic loci with equal effects.
In addition, there are 32 unlinked neutral loci with 256 possible
alleles subject to step-wise mutation \citep{oht73} at rate $10^{-3}$
per locus per generation. These loci are used to evaluate the levels
of genetic divergence within and between species that one would observe
if using microsatellite markers. Each of the parental species is initialized
with a specific allele at these loci to allow the origin of neutral
markers in the hybrids to be determined later.
\textbf{Life cycle.}\quad{}The life cycle consists of: (i) formation
of mating pairs, (ii) offspring production, (iii) viability selection
in eggs due to selection for local adaptation; (iv) selection in adults
due to predation and (v) dispersal. In \emph{Heliconius}, both selection
by predation and dispersal occur simultaneously. However, in our simulations
for computational purposes, these two processes were serialized. The
results presented below were obtained assuming selection before dispersal.
Limited simulations suggest that the alternative ordering leads to
similar results.
\textbf{Habitat preference.}\quad{}Following \citet{gav05a}, the
relative preference of an individual with habitat preference trait
$h$ for habitat $0$ (i.e. the open habitat) is given by a linear
function of $h$:
\begin{equation}
p=\epsilon+(1-2\epsilon)h,\label{eq:hab_pref}\end{equation}
where $\epsilon$ is a small positive parameter. The value of $p$
changes linearly from $\epsilon$ for $h=0$ to $1-\epsilon$ for
$h=1$. The value of $\epsilon$ can be interpreted as the probability
that an individual with a perfect preference for one habitat mistakenly
goes to the other habitat. The relative preference of this same individual
for habitat $1$ (i.e. the closed-canopy habitat) is $1-p$.
\textbf{Reproduction.}\quad{}Mating occurs between individuals in
the same {}``patch''. If a {}``patch'' has only one habitat type,
each male has an equal probability to encounter each female. If a
{}``patch'' has both habitat types, each male encounters a female
with a probability proportional to the similarity of their preferred
habitats. That is, a male and a female with habitat preference traits
$h_{m}$ and $h_{f}$ encounter or do not encounter each other with
probabilities $p_{m}p_{f}+(1-p_{m})(1-p_{f})$ and $p_{m}(1-p_{f})+p_{f}(1-p_{m})$,
respectively, where $p_{m}$ and $p_{f}$ are given by Equation~\ref{eq:hab_pref}
with an appropriate value of $h$. Note that the rate of hybridization
is an emergent character controlled by the current values of ecological
and preference traits.
Given an encounter, the actual mating occurs with a probability proportional
to the male's preference for the female's color pattern. The preference
of a male with preference traits $p_{r},p_{y}$ for a female with
color pattern $(i,j)$ is:
\begin{equation}
\psi(i,j\mid p_{r},p_{y})=\exp\left[\alpha_{r}(p_{r}-1/2)R_{ij}\right]\exp\left[\alpha_{y}(p_{y}-1/2)Y_{ij}\right],\label{psi}\end{equation}
where $R_{ij}$ and $Y_{ij}$ are the corresponding elements of the
color intensity matrices and $\alpha_{r}$ and $\alpha_{y}$ are positive
parameters measuring the strength of preference for red and yellow,
respectively. Large values of $\alpha$'s imply strong preferences;
small values imply weak preferences. This model represents a special
case of an open-ended mating preference (\citealt{lan82},\citealt[Chap.9]{gav04}).
Note that under our choice of the preference Function~\ref{psi},
each of the three species has the highest mating preference for its
own type \citep{nai01,mav06,mel09}.
Each mating results in a clutch of offspring the size of which is
drawn from a Poisson distribution with parameter $b$. We assume that
all adult females mate once. This assumption implies that any costs
of mate choice are absent, and that the effective population size
is increased relative to the actual number of adults \citep{gav05a}.
The clutch is laid on a host with a probability proportional to the
preference $p$ for host. Mutations occur at constant rate $\mu$
per gene per generation in all loci except for the sex-chromosome
$Z$ and the autosomal locus $C$ which do not mutate.
\textbf{Adaptation to host plants.}\quad{}In our model, adaptation
concerns the ability to grow on the host plant. The probability that
the clutch laid by a mother with an ecological trait $x$ survives
on a host is
\begin{equation}
w=\exp\left[-\frac{(x-\theta)^{2}}{2\sigma_{s}^{2}}\right],\end{equation}
where $\theta$ is the optimum phenotype, which is set to $0$ in
hosts of the open habitat, and to $1$ in hosts of the closed-canopy
habitat, and $\sigma_{s}$ is a positive parameter measuring the strength
of selection. This form of ecological selection introduces constraints
on local adaptation so that offspring cannot have high fitness on
both hosts. Smaller values of $\sigma_{s}$ imply stronger selection
for local adaptation and stronger selection against generalists (i.e.
individuals with ecological trait $x\sim0.5$) and specialists for
the alternative habitat. To account for the population size control
by the number of host plants available we assume that the number of
adults surviving selection for local adaptation is $K_{0}$ in habitat
0 and $K_{1}$ in habitat 1 where the carrying capacities $K_{0}$
and $K_{1}$ are parameters of the model.
\textbf{Selection by predation.}\quad{}Adults are subject to selection
by differential predation. We assume that predators (i.e. birds) stop
eating butterflies (which are distasteful) once they learn that a
particular color or a color pattern is associated with bad taste.
The efficiency of learning process depends both on the number of butterflies
eaten and the intensity of the stimuli (i.e. color and taste). We
posit that all butterflies are equally distasteful.
Two approaches are used here. The first assumes that birds learn to
avoid aposematic prey by separating the different color elements of
butterfly wing patterns. For a bird that has eaten $k$ butterflies,
we say its {}``learning scores'' for red and yellow colors are
\begin{equation}
Q_{r}=\sum_{j=1}^{k}R_{j},\ Q_{y}=\sum_{j=1}^{k}Y_{j},\end{equation}
where $R_{j}$ and $Y_{j}$ are the appropriate elements ({}``color
intensities'') of matrices $R$ and $Y$ for the $j$-th butterfly.
We assume that the bird stops eating the butterflies of a particular
color once the corresponding learning score reaches a positive threshold
$Q$. Small values of $Q$ imply that a small number of butterflies
will be eaten so that selection by predation is weak; large values
of $Q$ imply strong selection. The model assumes that more intense
colors (i.e. with larger values of $R$ and $Y$) are easier to remember.
This simple model is closely related to the ideas on the avoidance
of distasteful prey going back to \citet{mul79} (see \citealt{mal01}).
Note that once one color (red or yellow) is learned, selection on
the other color in mixed-color individuals disappears. We assume that
parameter $Q$ is the same for all patches. This assumption implies
that the density of predators is uniform across the whole system.
The second approach assumes that color patterns (i.e. combinations
of colors) are learned as a whole. By combining matrices $R$ and
$Y$, we get six unique patterns; a learning score $Q$ is introduced
for each of the six color patterns. The learning process happens as
before with each eaten butterfly contributing $1$ to the corresponding
learning score of a predator for the color pattern.
We note that the birds are not modeled explicitly but only implicitly
(via the learning thresholds) and that no evolution in birds was allowed.
\textbf{Dispersal.}\quad{}Individuals surviving predation disperse.
Butterflies are highly mobile and, thus, biologically realistic values
of the migration rates are large. We assume that each individual migrates
to one of the $8$ neighboring patches or stays in its native patch
with the same probability $1/9$. For the left-most and right-most
rows of patches which have 5 rather than 8 neighboring patches, a
proportion 3/9 of individuals are removed (reflecting their migration
outside the system).
\textbf{Initial conditions.}\quad{}The peripheral patches that contain
a single habitat type are initially populated with individuals perfectly
adapted to and with the highest preference for this habitat, and with
a genetic and phenotypic compositions of either \emph{H. melpomene}
(on the left side) or \emph{H. cydno} (on the right side). The patches
in the central area with both habitats were empty. 95\% of the individuals
are homozygous for all genes, but the remaining 5\% have an {}``intermediate''
sex chromosome (i.e. $Z_{2}$) to account for an incomplete segregation
between the two parental species, as explained above.
\textbf{Parameters.}\quad{}To analyze our model, we used individual-based
simulations, which we ran for $50,000$ generations, roughly corresponding
to $17,000$ years. We systematically varied four parameters: the
carrying capacity of habitat 1 ($K_{1}=200$ and $400$), clutch size
($b=4,8$, and $16$), the learning threshold ($Q=1,4$, and $16$)
and the strength of selection for local adaptation ($\sigma_{s}=0.1,0.25$,
and $0.4$). The following summarizes the parameters that did not
change: system size $96\times4$, the width of the central area 32
patches, the carrying capacity of habitat 0 $K_{0}=200$, mutation
rate $\mu=10^{-5}$ per generation, number of loci per trait $8$,
strength of color preference $\alpha_{r}=\alpha_{y}=4$, and probability
of habitat choosing error $\ge_{h}=0.001$. The reason for including
higher carrying capacity $K_{1}=400$ in habitat 1 is that some observations
suggest that \emph{H. cydno} has more hosts available than \emph{H.
melpomene} (Mavarez, pers.observation).
System-level data were saved every 100th generation; patch-level data
were saved every 1000th generation. The results below are based on
$20$ runs for each parameter combination. The model was implemented
in \emph{C}. The code is available upon request.
\section{Results}
After the start of each simulation run, individuals of both parental
species rapidly spread through the initially uninhabited central area
from the opposite sides and then start hybridizing, initially at low
frequencies. (The initial probability of hybridization was order $10^{-4}$
and the overall proportion of individuals with foreign neutral alleles
during the first 1000 generations was no larger than several percent.)
In the course of the simulations, we observed several different dynamic
regimes, and we will discuss them according to observed changes in
(i) coloration patterns, (ii) habitat use and local adaptation, and
(iii) reproductive compatibility. First, we discuss our results corresponding
to the case when birds learn different colors independently.
%
\begin{table}
\caption{\label{tab:results} {\small Numbers of different outcomes in the
dynamics of coloration patterns, ecological traits, and reproductive
compatibility (20 runs for each parameter combinations). {}``Red''
- allele $A$ fixed first, {}``yellow'' - allele $B$ fixed first.
{}``General'', {}``Tens'' and {}``Symp'' correspond to the formation
of the generalist ecotype, a tension zone separating two specialist
ecotypes, and two sympatrically distributed specialist ecotypes, respectively.
$Z_{2}$, $Z_{1}$ \& $Z_{3}$, and $Z_{2}$ \& $Z_{3}$ - signify
fixation of $Z_{2}$, to the loss of $Z_{2}$ and to the fixation
of $C_{2}$, respectively. The stars mark the combinations of parameters
discussed later in the text.}}
\begin{spacing}{0.8}
\hfill{} {\footnotesize }\begin{tabular}{|llll|cc|ccc|ccc|}
\hline
\multicolumn{4}{|c|}{{\footnotesize Parameters}} & \multicolumn{2}{|c|}{{\footnotesize Color}} & \multicolumn{3}{|c|}{{\footnotesize Ecotypes}} & \multicolumn{3}{|c|}{{\footnotesize $Z$ Chromosomes}}\tabularnewline
{\footnotesize $K_{1}$ } & {\footnotesize $b$ } & {\footnotesize $Q$ } & {\footnotesize $\sigma_{s}$ } & {\footnotesize Red } & {\footnotesize Yellow } & {\footnotesize General } & {\footnotesize Tens } & {\footnotesize Symp } & {\footnotesize $Z_{2}$ fixed } & {\footnotesize $Z_{2}$ lost } & {\footnotesize $C_{2}$ fixed }\tabularnewline
\hline
{\footnotesize 200 } & {\footnotesize 4.0 } & {\footnotesize 1 } & {\footnotesize 0.10 } & {\footnotesize 13 } & {\footnotesize 7 } & {\footnotesize 0 } & {\footnotesize 13 } & {\footnotesize 7 } & {\footnotesize 14 } & {\footnotesize 6 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25 } & {\footnotesize 8 } & {\footnotesize 12 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 16 } & {\footnotesize 4 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40 } & {\footnotesize 5 } & {\footnotesize 15 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 17 } & {\footnotesize 3 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & {\footnotesize 4 } & {\footnotesize 0.10 } & {\footnotesize 10 } & {\footnotesize 10 } & {\footnotesize 0 } & {\footnotesize 17 } & {\footnotesize 3 } & {\footnotesize 18 } & {\footnotesize 2 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25 } & {\footnotesize 14 } & {\footnotesize 6 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 15 } & {\footnotesize 4 } & {\footnotesize 1 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40 } & {\footnotesize 9 } & {\footnotesize 11 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 19 } & {\footnotesize 0 } & {\footnotesize 1 }\tabularnewline
\noalign{\vskip-2pt}
& & {\footnotesize 16 } & {\footnotesize 0.10 } & {\footnotesize 11 } & {\footnotesize 9 } & {\footnotesize 0 } & {\footnotesize 18 } & {\footnotesize 2 } & {\footnotesize 17 } & {\footnotesize 2 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25 } & {\footnotesize 11 } & {\footnotesize 9 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 18 } & {\footnotesize 1 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40 } & {\footnotesize 15 } & {\footnotesize 5 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 18 } & {\footnotesize 1 } & {\footnotesize 1 }\tabularnewline
\noalign{\vskip-2pt}
& {\footnotesize 8.0 } & {\footnotesize 1 } & {\footnotesize 0.10$^{*}$ } & {\footnotesize 9 } & {\footnotesize 11 } & {\footnotesize 0 } & {\footnotesize 2 } & {\footnotesize 18 } & {\footnotesize 13 } & {\footnotesize 7 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25$^{*}$ } & {\footnotesize 4 } & {\footnotesize 16 } & {\footnotesize 0 } & {\footnotesize 17 } & {\footnotesize 3 } & {\footnotesize 17 } & {\footnotesize 3 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40$^{*}$ } & {\footnotesize 9 } & {\footnotesize 11 } & {\footnotesize 18 } & {\footnotesize 2 } & {\footnotesize 0 } & {\footnotesize 18 } & {\footnotesize 2 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & {\footnotesize 4 } & {\footnotesize 0.10 } & {\footnotesize 14 } & {\footnotesize 6 } & {\footnotesize 0 } & {\footnotesize 4 } & {\footnotesize 16 } & {\footnotesize 18 } & {\footnotesize 2 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25 } & {\footnotesize 8 } & {\footnotesize 12 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40 } & {\footnotesize 12 } & {\footnotesize 8 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & {\footnotesize 16 } & {\footnotesize 0.10 } & {\footnotesize 16 } & {\footnotesize 4 } & {\footnotesize 0 } & {\footnotesize 2 } & {\footnotesize 18 } & {\footnotesize 18 } & {\footnotesize 2 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25 } & {\footnotesize 12 } & {\footnotesize 8 } & {\footnotesize 0 } & {\footnotesize 19 } & {\footnotesize 1 } & {\footnotesize 18 } & {\footnotesize 1 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40 } & {\footnotesize 14 } & {\footnotesize 6 } & {\footnotesize 19 } & {\footnotesize 1 } & {\footnotesize 0 } & {\footnotesize 19 } & {\footnotesize 1 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& {\footnotesize 16.0 } & {\footnotesize 1 } & {\footnotesize 0.10 } & {\footnotesize 11 } & {\footnotesize 9 } & {\footnotesize 0 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 19 } & {\footnotesize 1 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25 } & {\footnotesize 9 } & {\footnotesize 11 } & {\footnotesize 5 } & {\footnotesize 13 } & {\footnotesize 2 } & {\footnotesize 18 } & {\footnotesize 2 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40 } & {\footnotesize 6 } & {\footnotesize 14 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 } & {\footnotesize 16 } & {\footnotesize 4 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & {\footnotesize 4 } & {\footnotesize 0.10 } & {\footnotesize 11 } & {\footnotesize 9 } & {\footnotesize 0 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 17 } & {\footnotesize 3 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25 } & {\footnotesize 8 } & {\footnotesize 12 } & {\footnotesize 9 } & {\footnotesize 11 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40 } & {\footnotesize 10 } & {\footnotesize 10 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & {\footnotesize 16 } & {\footnotesize 0.10 } & {\footnotesize 13 } & {\footnotesize 7 } & {\footnotesize 0 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 19 } & {\footnotesize 1 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25 } & {\footnotesize 13 } & {\footnotesize 7 } & {\footnotesize 7 } & {\footnotesize 13 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40 } & {\footnotesize 9 } & {\footnotesize 11 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\hline
{\footnotesize 400 } & {\footnotesize 4.0 } & {\footnotesize 1 } & {\footnotesize 0.10 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 8 } & {\footnotesize 3 } & {\footnotesize 9 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 11 } & {\footnotesize 3 } & {\footnotesize 6 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 16 } & {\footnotesize 1 } & {\footnotesize 3 }\tabularnewline
\noalign{\vskip-2pt}
& & {\footnotesize 4 } & {\footnotesize 0.10 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 5 } & {\footnotesize 1 } & {\footnotesize 14 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 13 } & {\footnotesize 1 } & {\footnotesize 6 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & {\footnotesize 16 } & {\footnotesize 0.10 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 2 } & {\footnotesize 4 } & {\footnotesize 14 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 12 } & {\footnotesize 0 } & {\footnotesize 8 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& {\footnotesize 8.0 } & {\footnotesize 1 } & {\footnotesize 0.10 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 17 } & {\footnotesize 0 } & {\footnotesize 3 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25$^{*}$ } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 19 } & {\footnotesize 0 } & {\footnotesize 1 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 19 } & {\footnotesize 1 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & {\footnotesize 4 } & {\footnotesize 0.10 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 18 } & {\footnotesize 0 } & {\footnotesize 2 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & {\footnotesize 16 } & {\footnotesize 0.10 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 19 } & {\footnotesize 0 } & {\footnotesize 1 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40 } & {\footnotesize 1 } & {\footnotesize 19 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& {\footnotesize 16.0 } & {\footnotesize 1 } & {\footnotesize 0.10 } & {\footnotesize 3 } & {\footnotesize 17 } & {\footnotesize 0 } & {\footnotesize 16 } & {\footnotesize 4 } & {\footnotesize 17 } & {\footnotesize 3 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 16 } & {\footnotesize 4 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & {\footnotesize 4 } & {\footnotesize 0.10 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 16 } & {\footnotesize 4 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40 } & {\footnotesize 1 } & {\footnotesize 19 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & {\footnotesize 16 } & {\footnotesize 0.10 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 17 } & {\footnotesize 3 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.25 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\noalign{\vskip-2pt}
& & & {\footnotesize 0.40 } & {\footnotesize 1 } & {\footnotesize 19 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 20 } & {\footnotesize 0 } & {\footnotesize 0 }\tabularnewline
\hline
\end{tabular}\hfill{}\end{spacing}
\end{table}
%
\begin{figure}[t]
\hfill{}\includegraphics[width=2in]{Figures/histo0}\includegraphics[width=2in]{Figures/histo5000}\includegraphics[width=2in]{Figures/histo10000}\hfill{}
\hfill{}\includegraphics[width=2in]{Figures/histo15000}\includegraphics[width=2in]{Figures/histo20000}\includegraphics[width=2in]{Figures/histo25000}\hfill{}
\caption[~~~The dynamics of the frequencies of \textbf{A} (red) and \textbf{B}
(yellow) alleles across the whole system]{ The dynamics of the frequencies of \textbf{A} (red) and \textbf{B}
(yellow) alleles across the whole system. The figure includes all
runs for the first 25,000 generations.}
\label{colors}
\end{figure}
\textbf{The dynamics of coloration patterns.}\quad{}In the simulations
there was a strong tendency for fixation of both red (A) and yellow
(B) alleles across the whole system so that a \emph{Heurippa}-like
color pattern eventually dominated (see Figure~\ref{colors}). These
dynamics are explained by two factors. First, natural selection by
predation favors the establishment of aposematic coloration and, thus,
increases the frequencies of the warning coloration alleles. Second,
since parental species already have mating preference alleles for
the corresponding colors, hybrids acquire them as well. As a result,
sexual selection favors the fixation of the {}``red'' and {}``yellow''
alleles via a Fishering runaway process \citep{lan81,and94}. In fact,
the color preference traits $p_{r}$ and $p_{y}$ approach the maximum
possible value of one across the whole system much faster than either
of the color alleles approach fixation: on average, in 3,600 generations.
In contrast, the fixation of a first color allele (red or yellow)
takes place on average in 6,300 generations. The second color locus
approaches fixation much later - on average about 16,100 generations
later. The reason for this discrepancy is that while the initial evolution
of the color loci is driven both by natural and sexual selection,
once one color locus is close to fixation, natural selection by predation
is not acting on the second locus anymore and its evolution is mostly
driven by mating preferences (and random genetic drift).
Whichever color locus is fixed first controls for which locus the
polymorphism will be transiently maintained and, thus, can be observed
in natural populations (see Table~\ref{tab:results} for the frequencies
at which the red and yellow alleles are first to be fixed for different
parameter values). If the carrying capacities of both hosts are equal
(i.e. $K_{0}=K_{1}=200$), fixation of the red allele (and, consequently,
long lasting polymorphism in the yellow locus) is favored by an intermediate
clutch size ($b=8.0$). It is also favored by strong selection for
local adaptation (i.e. low values of $\sigma_{s}$) coupled with extreme
values of the learning threshold $Q$ ($Q=1$ and $16$). If the carrying
capacity of habitat 1 is twice that of habitat 0 (i.e. $K_{1}=400$),
the yellow allele is almost always fixed first while the red locus
remains polymorphic. This is expected because initially the system
as a whole has twice as many $B$ alleles as $A$ alleles.
\textbf{The dynamics of ecological traits.}\quad{}Populations inhabiting
the peripheral areas with only one host type present always remain
close to the original values of the ecological traits (and, correspondingly,
remain locally adapted) independently of the strength of selection
for local adaptation $\sigma_{s}$. The situation, however, is different
in the central area, where both types of hosts are available for butterflies
to mate and oviposit. Three different long-term outcomes are observed.
First, most or all individuals in the central area have intermediate
values of the ecological trait close to $0.5$ (Figure~\ref{ecology}a).
These are the generalists utilizing both host species. At the boundaries
of the central and peripheral areas, two tension zones are formed
separating the populations of generalists in the central area and
specialists in the peripheral areas. {[}A tension zone as a hybrid
zone maintained by a balance of dispersal of parental forms and selection
against hybrids.{]} Second, a single tension zone is formed in the
central area which then moves and eventually stabilizes at a boundary
of the central area. This tension zone separates two specialist populations
and largely prevents the utilization of one host in the central area
(Figure~\ref{ecology}b). Finally, two specialist populations (each
adapted to and utilizing its own host) coexist sympatrically through
the central area (Figure~\ref{ecology}c). Note that in this case
the traits values in the central area are shifted away from the optimum
values because of the gene flow between the specialists due to hybridization.
In all three cases, the habitat preference trait $h$ closely matches
the ecological trait across the whole system.
%
\begin{figure}[t]
\hfill{}\subfloat[]{
\includegraphics[width=2in]{Figures/habitat1BW}}\subfloat[]{
\includegraphics[width=2in]{Figures/habitat2BW}}\subfloat[]{
\includegraphics[width=2in]{Figures/habitat3BW}}\hfill{}
\caption[~~~Three different outcomes in the dynamics of the average ecological
character $x$ per habitat across space]{\label{ecology} Three different outcomes in the dynamics of the
average ecological character $x$ per habitat across space. (a)~Central
area is inhabited by a population of generalists. (b)~Two specialists
separated by a tension zone at a border of the central area. (c)~Two
specialists coexisting throughout the central area. The gray area
marks the central patches. Parameter values used to produce the figures:
(a)~$K_{1}=200,b=16,Q=1,\sigma_{s}=0.4$, (b)~$K_{1}=200,b=8,Q=4,\sigma_{s}=0.25$,
(c)~$K_{1}=200,b=8,Q=4,\sigma_{s}=0.1$. The two lines in the central
region show average trait values computed over two hosts.}
\end{figure}
Consider first the situation when the carrying capacities of both
hosts are equal (i.e. $K_{0}=K_{1}=200$). In this case, the most
common outcome is the formation of a tension zone separating two specialists
(see Table \ref{tab:results}, Figure~\ref{ecology}b). Generalists
(as in Figure~\ref{ecology}a) are more likely to emerge if the birth
rate $b$ and the coefficient $\sigma_{s}$ are large so that selection
for local adaptation is relatively weak. Sympatric coexistence of
two specialists in the central area (as in Figure~\ref{ecology}c)
was observed only few times and required selection for local adaptation
to be the strongest (that is, $\sigma_{s}=0.1$). If the carrying
capacity of habitat 1 is twice that of habitat 0 (i.e. $K_{1}=400$),
generalists were never observed. Two sympatric specialists coexisting
in the central area were observed only a couple of times, leaving
the formation of a tension zone separating two specialists (Figure~\ref{ecology}b)
as the most common outcome.
\textbf{The dynamics of reproductive compatibility.}\quad{}Recall
that the simulations started with two parental species almost fixed
for the reproductively isolated combinations of sex chromosomes and
autosomal genes $Z_{1}C_{1}C_{1}$ and $Z_{3}C_{2}C_{2}$, respectively,
with only a small proportion of the \emph{H. heurippa}-type chromosome
$Z_{2}$ present in both parental species. At the end of the simulations,
four situations were observed. First, $Z_{2}$ can be lost, with the
reproductively isolated ancestral combinations $Z_{1}C_{1}C_{1}$
and $Z_{3}C_{2}C_{2}$ being both present in the system. Second, \emph{H.
heurippa}-type $Z_{2}$ can be fixed across the whole system, while
the autosomal locus $C$ can exhibit neutral polymorphism with both
$C_{1}$ and $C_{2}$ alleles present. Third, the \emph{H. melpomene}-type
autosomal allele $C_{1}$ can be fixed while sex chromosomes $Z_{1}$
and $Z_{2}$ are both maintained in a neutral fashion. Finally, the
\emph{H. cydno}-type autosomal allele $C_{2}$ can be fixed while
sex chromosomes $Z_{3}$ and $Z_{2}$ are both maintained in a neutral
fashion. In the last three cases, all genotypes present in the system
are reproductively compatible.
If the carrying capacities of both hosts are equal (i.e. $K_{0}=K_{1}=200$),
the most common outcome is the fixation of \emph{H. heurippa}-type
sex chromosome $Z_{2}$ in the whole system (about 90\% of cases).
The loss of $Z_{2}$ which was the second most common outcome (almost
9\% of cases) was promoted by smaller birth rate $b$ and stronger
selection for local adaptation (i.e. smaller $\sigma_{s}$). After
the loss of the $Z_{2}$ chromosome, a future hybrid speciation event
is unlikely. The fixation of \emph{H. melpomene}-type allele $C_{1}$
or \emph{H. cydno}-type allele $C_{2}$ with a neutral polymorphism
in the $Z$ chromosome were observed only a handful of times.
If the carrying capacity of habitat 1 is twice that of habitat 0 (i.e.
$K_{1}=400$), the most common outcome again was the fixation of $Z_{2}$.
However now with a small birth rate (i.e. $b=4$), the fixation of
\emph{H. cydno}-type allele $C_{2}$ with the neutral coexistence
of $Z_{3}$ and $Z_{2}$ was frequently observed for strong and intermediate
strengths of selection for local adaptation (i.e. $\sigma_{s}=0.1$
and $0.25$). This is what is expected because of higher initial number
of \emph{H. cydno}-type alleles in the system. The loss of $Z_{2}$
was observed few times (mostly with $b=4$) while the fixation of
$C_{1}$ with the neutral coexistence of $Z_{1}$ and $Z_{2}$ was
never observed.
\textbf{Joint dynamics.}\quad{}The hybrid aposematic coloration (i.e.
the presence of both red and yellow bands) readily spreads across
the whole system in all cases considered. This spread is preceded
by the establishment of strong mating preferences for both red and
yellow colors. The spread of hybrid coloration is most often accompanied
by the fixation of the hybrid chromosome $Z_{2}$ in the whole system
so that all individuals are reproductively compatible. With regard
to the ecological and habitat preference characters, the two ecotypes
initially present in the parental species (which have high preference
for and local adaptation on hosts in habitat 0 and 1, respectively)
are preserved and can either coexist thorough the central area or
one of them occupied most of this area. It is also possible that an
additional generalist ecotype utilizing both habitats can emerge and
spreads across the central area excluding both specialists.
The loss of $Z_{2}$ results in the presence of two genetically incompatible
sets of populations with \emph{H. melpomene}-type genotype $Z_{1}C_{1}C_{1}$
and \emph{H. cydno}-type genotype $Z_{3}C_{2}C_{2}$ both having \emph{H.
heurippa}-type coloration and mating preferences. Most commonly these
populations are separated by a tension zone located at a border of
the central area but in a few cases they are sympatric in this area.
In few other cases, one of the two sets of populations forms a generalist
ecotype that spreads across the central area excluding populations
of the other type and the specialist ecotype of its own type.
The fixation of \emph{H. cydno}-type autosome $C_{2}$ occurring simultaneously
with the neutral coexistence of $Z_{2}$ and $Z_{3}$ was accompanied
by the separation of the two specialist ecotypes by a tension zone
located at a border of the central area.
Figure~\ref{fig:timescale} shows the typical timescales in our simulations.
If \emph{H. heurippa}-type chromosome $Z_{2}$ is destined to be lost,
this happens relatively fast, typically during the second thousand
of generations. Independently of this, within a few thousand generations
after the start of the simulation both mating preference traits evolve
to extreme values and the ecological and the habitat preference traits
approach their equilibrium distributions. Few thousand generations
later, one of the color pattern alleles approaches fixation and the
\emph{H. heurippa}-type sex chromosome $Z_{2}$ achieves high frequency
(if $Z_{2}$ is destined to be fixed). The other color allele locus
approaches fixation far more slowly (data not shown). The growth in
the frequencies of the color alleles and $Z_{2}$ chromosome start
at the same border of the central area and then spreads throughout
the whole system.
Technical comment: We note that here the time for a color allele to
be {}``fixed'' means the time for this allele to reach a frequency
of 95\% in the whole system. {[}We use 95\% rather than 100\% because
of the recurrent mutation.{]} The time for a particular sex chromosome
to be {}``fixed'' or {}``lost'' means the time for the chromosome
to reach a frequency of 100\% or 0\%. To make conclusions about ecological
traits we use a more convoluted procedure. First, we identify all
patches in the central area in which each host plant contributes at
least 20\% of the patch's population. We say that a {}``generalist''
is formed if among these patches there is a band of patches of at
least 6 patches wide in which the trait values are close to $1/2$
(specifically, the average ecological trait of the individuals born
in environment type 0 is above 0.4 and for those born in type 1 is
below 0.6). We say that there is sympatric coexistence of two ecomorphs
if there is a band of patches at least 6 patches wide in which the
trait values are away from $1/2$ (specifically, the average ecological
trait of the individuals born on host 0 is below 0.4 and for those
born on host 1 is above 0.6). We say there is a {}``tension zone''
otherwise. That is, either the width of the band of patches with two
populations is 5 patches or less, or there are demes with high average
ecological character (on one side of the tension zone) and demes with
low average value (on the other side).
%
\begin{figure}[t]
\hfill{}\subfloat[]{
\includegraphics[width=3in]{Figures/timescale_200}}\subfloat[]{
\includegraphics[width=3in]{Figures/timescale_400}} \hfill{}
\caption[~~~Tukey plots of the timescales associated with different outcomes
marked on the $x$-axis]{ \label{fig:timescale} Tukey plots of the timescales associated
with different outcomes marked on the $x$-axis. The y-axis represents
time in generations on a logarithmic scale. The vertical lines extend
from minimum to maximum observations, the middle dashed lines depict
averages, and the boxes extend from lower to upper quartiles. (a)~$K_{1}=200$.
(b)~$K_{1}=400$.}
\end{figure}
\textbf{The dynamics of neutral markers.}\quad{}The above dynamics
were accompanied by the introgression in neutral markers. Figure~\ref{fig:markers}
shows the frequency of neutral markers of \emph{H. melpomene} origin
across the whole system estimated by randomly sampling 16 individuals
from each of the 96 patches in the first row of patches. We used four
parameter combinations which showed a diversity of outcomes in the
dynamics of coloration, habitat usage, and reproductive compatibility
(see Table~ \ref{tab:results} where the combinations of parameters
used are marked with a star). One can see that at the onset of hybridization
the frequency of \emph{H. melpomene} markers experiences a rapid shift
away from its initial value of $0.5$. This shift is apparently explained
by strong sexual selection operating during a relatively short time
interval when mating preference for one color has already reached
the highest value while the preference for the second color is still
intermediate. During this time interval, individuals with the maximum
expression of the second color have mating advantage which results
in increasing the frequency of their neutral markers in the population.
The explanation just given works for the data shown in Figures~\ref{fig:markers}a,b
and d but not for those in Figures~\ref{fig:markers}c. This suggests
that some additional factors are also in place. Note that the downward
shift in the frequency of neutral markers of \emph{H. melpomene} origin
observed in Figures~\ref{fig:markers}d is explained by the fact
that for this parameter combination initially there are twice as many
\emph{H. cydno} individuals as \emph{H. melpomene }individuals. After
both mating preference traits have reached the maximum value of one,
the dynamics of neutral markers appear to be driven mostly by random
genetic drift.
%
\begin{figure}[p]
\hfill{}\subfloat[]{
\includegraphics[width=3in]{Figures/plot1}}\subfloat[]{
\includegraphics[width=3in]{Figures/plot2}}\hfill{}
\hfill{}\subfloat[]{
\includegraphics[width=3in]{Figures/plot3}}\subfloat[]{
\includegraphics[width=3in]{Figures/plot4}}\hfill{}
\caption[~~~The dynamics of the frequency of neutral markers of \emph{H.
melpomene} origin for four parameter combinations]{\label{fig:markers} The dynamics of the frequency of neutral markers
of \emph{H. melpomene} origin for four parameter combinations (marked
with a star in Table~1). In all figures, $Q=1$ and $b=8.0$. (a)~$K_{1}=200,\sigma_{s}=0.1$
(b)~$K_{1}=200,\sigma_{s}=0.25$ (c)~$K_{1}=200,\sigma_{s}=0.4$
(d)~$K_{1}=400,\sigma_{s}=0.25$. 20 runs for each combination of
parameters. Solid lines: runs in which the preference for yellow color
trait $p_{y}$ reached $0.99$ first. Dashed lines: runs in which
the preference for red color trait $p_{r}$ reached $0.99$ first.}
\end{figure}
\afterpage{\clearpage}
Figure~\ref{neut}a gives more details on the dynamics of neutral
markers using a numerical run with parameter values as in Figure~\ref{fig:markers}a.
Shown is the whole distribution of the frequency of the neutral markers
of \emph{H. melpomene} origin (rather than just its average value
as Figure~\ref{fig:markers}). The distribution becomes very wide
soon after the onset of hybridization. Figure~\ref{neut}b shows
the corresponding dynamics of the frequencies of the neutral markers
of \emph{H. melpomene} origin estimated in three groups of four demes
each of the first row: four leftmost demes, four centrally located
demes and four rightmost demes. This figure shows that it takes a
long time for neutral markers to spread across the central area and
then to reach the most peripheral patches.
%
\begin{figure}[t]
\hfill{}\subfloat[]{
\includegraphics[width=2.9in]{Figures/contour}}\subfloat[]{
\includegraphics[width=2.9in]{Figures/3regions}}\hfill{}
\caption[~~~The dynamics of neutral markers of \emph{H. melpomene} origin
in a run with $K_{1}=200,Q=1,b=8.0$ and $\sigma_{s}=0.1$]{\label{neut} The dynamics of neutral markers of \emph{H. melpomene}
origin in a run with $K_{1}=200,Q=1,b=8.0$ and $\sigma_{s}=0.1$.
(a)~The distribution of the frequency of neutral markers of \emph{H.
melpomene} origin over time. The intensity of the black color is proportional
to the number of individuals with the corresponding frequency (the
scale is given by the color bar on the right). (b)~The average frequencies
of neutral markers of \emph{H. melpomene} origin in four leftmost,
four centrally located, and four rightmost demes.}
\end{figure}
\textbf{Some transient dynamics.}\quad{}In no case did we observe
a stable coexistence of the three color patterns of parental and hybrid
species. Was there some kind of transient coexistence which did not
persist to the end of our simulations? To answer this question we
identified time intervals and demes in which the two parental and
a hybrid phenotypes where present at frequencies higher than 10\%
(chosen to be larger than the maximum observed proportion of natural
hybrids of 8\%). Across all runs, this transient period starts on
average around generation 2,200 (std. dev. 1,200) and no later than
generation 5,000. It lasts on average 900 generations (std. dev. 800),
with a maximum of 5,000 generations. An analysis of variances shows
that both the starting generation and duration of the transient period
decrease with increasing the strength of predation (i.e. decreasing
the value of $Q$). The duration of the transient period is further
affected by the clutch size; extreme values of $b$ (i.e. $4$ and
$16$) promote longer transients.
\textbf{Color patterns learned as a whole}.\quad{}In this case, novel
combinations of colors are under strong negative selection by predation.
As a result, hybrids never achieve high frequencies, $Z_{2}$ chromosome
disappears from the system, and both parental species are preserved
and coexist sympatrically in the central area (see Figure \ref{ecology}c)
for all combinations of parameters considered (data not shown).
\section{Discussion}
In this paper we used a spatial individual-based multilocus model
to study a proposed case of homoploid hybrid speciation involving
two sister species, \emph{Heliconius cydno} and \emph{H. melpomene},
and a putative hybrid species \emph{H. heurippa} \citep{mav06}. Assessing
the mechanisms potentially driving hybrid speciation and the plausibility
of its different scenarios can help us better understand both the
origins and evolution of \emph{Heliconius} butterflies and the challenges
faced by empiricists studying homoploid hybrid speciation in animals.
Our results show that if the birds learn different colors independently,
then \emph{H. heurippa}-type coloration pattern and the mating preferences
for this pattern as well as \emph{H. heurippa}-type sex chromosome
can readily spread across the whole system within ten thousand generations
after the initial contact between the parental species. The resulting
hybrid species retains the two ecotypes of the parental species using
two different hosts and can also develop an additional generalist
ecotype utilizing both hosts. The emerging hybrid species has neutral
markers from both parental species at substantial frequencies. It
is also possible that the two parental species converge on the same
hybrid coloration and mating preferences and undergo substantial exchange
in neutral markers while retaining their original sex chromosomes
and ecotypes. This outcome can be interpreted as the evolution of
Mullerian mimicry as a result of hybridization. However, empirical
data provide little support for this alternative outcome given that
most \emph{Heliconius} sister species pairs (i.e. able to hybridize)
tend to belong to different mimicry rings. In contrast, species in
the same mimicry ring are usually very distantly related and unable
to hybridize \citep{bel00}.
The spread of hybrid coloration is largely driven by selection for
warning coloration by predation. It is also helped by the pre-existing
(in the parental species) mating preferences for bright colors which
pass to the hybrids \citep{mel09} and then coevolve with coloration
patterns in a process analogous to Fisherian runaway. In the \emph{Heliconius
}butterflies, coloration plays the role of a {}``magic trait'' \citep{gav04,kro06,jig08,jig08a,sal08,pue07}
that controls both viability and mating. The existence of magic traits
greatly simplifies non-allopatric speciation in general \citep{gav04,pue07}
and in butterflies in particular \citep{mal07,jig08,jig08a,sal08}.
The evolution of ecotypes is driven by selection for local adaptation.
If this selection is weak, a generalist ecotype utilizing both host
can emerge in the central area and expel the specialist ecotypes from
this area (c.f., \citealt{gav07a,gav07b}). If selection is strong,
the gene flow from one specialist to another is small and does not
prevent the colonization of the whole central area by both ecotypes.
For intermediate strengths of selection, the formation of the generalist
ecomorph is thwarted and simultaneously the gene flow between the
two specialist ecomorphs is strong enough to prevent their coexistence
in the patches located in the central area. As a result, one specialist
ecomorphs largely excludes the other one from the central area. This
happens as the tension zone moves stochastically across the central
area until it gets trapped at its boundary where selection regime
changes \citep{bar83,bar99}. The explanation of why \emph{H. heurippa}-type
sex chromosome $Z_{2}$ spreads through the system is a bit more convoluted.
After the mating preferences for both red and yellow colors have spread
through the system, mating becomes random with respect to sex chromosomes.
In this case, while individuals with \emph{H. melpomene}-type chromosome
$Z_{1}$ and \emph{H. cydno}-type chromosome $Z_{3}$ are mutually
incompatible (in the sense that female hybrid offspring is infertile),
those with $Z_{2}$ are partially compatible with both parental types.
As a result, individuals with $Z_{2}$ chromosome have on average
more fertile offspring than those with either $Z_{1}$ or $Z_{3}$.
This ultimately drives $Z_{2}$ chromosome towards fixation.
In contrast, if different color patterns are learned as a whole, then
selection by predation prevents any significant accumulation of hybrids
in the system and hybrid speciation does not occur.
Recently, \citet{sal08} and \citet{jig08a} have contrasted two modes
of homoploid hybrid speciation which they called \emph{hybrid trait
speciation} and \emph{mosaic genome hybrid speciation}. In the former
mode, speciation occurs through the establishment via hybridization
of a novel adaptive trait. The novel trait must also confer a degree
of reproductive isolation from the parental lineages. In the later
mode, speciation involves stabilization of a hybrid genome which initially
contains a large number of intrinsically incompatible genes. The two
modes, thus, mostly differ in the type of selection involved (that
is, ecological vs. intrinsic incompatibility). \citet{jig08a} also
expected them to differ in the proportion of the genome involved with
hybrid trait speciation resulting in introgression of a handful of
genes and in mosaic genome hybrid speciation resulting in massive
introgression (as in the case of sunflowers, e.g. \citealt{ung98}).
According to \citet{sal08} and \citet{jig08a} who analyzed genetic
divergence between \emph{H. melpomene}, \emph{H. cydno} and \emph{H.
heurippa} in several loci, this case fits the hybrid trait speciation
mode. Our data do not support the argument of \citet{sal08} and \citet{jig08a}.
In our model, hybrid speciation is indeed mostly driven by evolutionary
processes involving a single trait - coloration. However, the extent
of hybridization in neutral loci is very extensive and is more compatible
with that expected under mosaic genome hybrid speciation. This suggests
that discriminating between hybrid trait speciation and mosaic genome
hybrid speciation may not be easy.
A question of the speed of ecological speciation is of great importance
and interest \citep{hen07,gav07a,gav07b,thi08} with most models suggesting
that speciation (i.e. evolution of significant reproductive isolation)
could happen on the time scale of a few hundred to a few thousand
generations. The data in \citet{mav06} have sometimes been interpreted
as suggesting that hybrid speciation could happen within a few generations.
Our model does not support this interpretation. Indeed in our simulations,
the typical time-scales of the evolution of Mullerian mimicry and
of hybrid speciation are on the order of several thousand generations.
In our simulations, the initial rate of hybridization is very low.
If it were higher, the hybrid species would evolve faster. The rate
of hybrid speciation is expected to increase if the hybrids utilize
a novel habitat \citep{nol05,schw07} that is free of the parental
forms.
In our simulations we never observed the coexistence of the two parental
and the hybrid species. The same mechanisms that facilitated the emergence
and establishment of the hybrid species also lead to the spread of
the hybrid species traits across the whole system. As a result, the
parental combinations of traits disappear (colors always, $Z$ chromosomes
most of the time). The same outcome is expected after a secondary
contact between the hybrid and a parental species. Therefore, our
model suggests that the coexistence of the three \emph{Heliconius}
species in South America may be transient rather than stable. An alternative,
and likely more reasonable, explanation is that our model lacks certain
components crucial for the coexistence of the three species. One possibility
is spatial and temporal heterogeneity of the biotic environment which
can result in a partial fragmentation of the species ranges (e.g.
if the density of hosts varies dramatically in space) and/or spatial
variation of selection (both for adaptation to hosts and by predation).
This heterogeneity will reduce the opportunity for gene flow and might
ensure coexistence. Alternatively, selection by predation can act
via some other mechanism that would increase geographic variation
in the system. For example, it is possible that parental colorations
are protected by selection if both parental species have some additional
Mullerian mimics, which is indeed the case in Heliconiinae butterflies.
Evolutionary dynamics of Mullerian mimicry can be complex especially
if different members of a Mullerian ring differ in the degree of unpalatability
\citep{gav98d}. This may further complicate the dynamics of hybrid
speciation and make predictions difficult. Other factors that may
affect the dynamics are selection for local adaptation and the genetic
architecture of the traits considered. In the model we made the simplest
assumptions of additive genetics and unlinked genes. More complex
schemes can of course alter the dynamics.
Our model also predicts the appearance of \emph{H. heurippa}-type
individuals utilizing \emph{H. melpomene}-type hosts and of a generalist
ecotype utilizing both hosts. Apparently, neither of the these two
outcomes has been observed in nature. While evolution of the generalist
is predicted to happen only under some conditions (specifically, under
weak selection for local adaptation), some utilization of the open
habitat initially utilized by \emph{H. melpomene} occurs in the model
always. It is well known that although tension zones represent a barrier
to the neutral genes (e.g. \citealt{ben85,bar86,gav97a,gav98d}),
advantageous genes will pass through the tension zones easily. What
factors would prevent the spread of \emph{H. heurippa} coloration
across \emph{H. melpomene}-like ecotype once hybridization starts
is unclear but one possibility is that the signaling quality of aposematic
color patterns varies with the habitat type or the environmental background
\citep{swe03}. The \emph{H. heurippa} colour pattern might therefore
be less fit to the open habitat areas typically occupied by \emph{H.
melpomene}.
Unfortunately, the lack of data does not yet allow one to come up
with a mathematical model that would have a better explanatory power.
More and better data are needed to that end. From the theoretical
perspective, at this stage the most crucial would be to understand
better how selection by predation operates in the system, how birds
learn the colors, and to what extent they distinguish different components
of color patterns. Of course any data on genetics of the traits considered
including those controlling reproductive compatibility would greatly
increase the realism of the model. The same effect will be achieved
by more precise information of the spatial distributions and densities
of the parental species and their co-mimics, and on the degree of
temporal stability of \emph{H. heurippa} geographic distribution.
Without such data it will remain difficult to make definite conclusions
on the origins of \emph{H. heurippa}.
However, one can attempt to identify possible scenarios from both
empirical data and mathematical models. \emph{H. heurippa} is well
differentiated from both \emph{H. melpomene} and \emph{H. cydno} in
some microsatellite loci (Mavarez, unpubl.data). In addition, the
current range of \emph{H. heurippa} is well outside of the area of
sympatry of \emph{H. melpomene} and \emph{H. cydno}. If one accepts
the hybrid origin of \emph{H. heurippa}, these two observations suggest
both long-lasting geographic isolation of \emph{H. heurippa} from
its parental species and significant changes in the species ranges
since speciation. Overall, our model supports the possibility of hybrid
origin of \emph{H. heurippa}. The most plausible scenario would include
hybridization between \emph{H. melpomene} and \emph{H. cydno} in an
area geographically isolated from the rest of both parental species
with subsequent long-lasting geographic isolation of the new hybrid
species, followed by changes in the species ranges, the secondary
contact, and the disappearance of \emph{H. melpomene}-like ecomorph
in the hybrid species. Our model does not rule out completely an alternative
scenario of the appearance of the red allele \textbf{A} in a geographically
isolated \emph{H. cydno} population by mutation, with the subsequent
fixation of the \emph{H. heurippa}-type color. However, this scenario
appears less likely at least for two reasons. First, it would require
a very rare mutational event that would bring the red allele. Second,
it would require an additional sequence of unlikely events that would
fix this red allele by drift in the presence of positive frequency-dependent
selection due to predation acting to reduce the frequency of this
initially rare allele \citep{mal89}. Although such stochastic peak
shifts can occur, their probability is very small \citep{bar84,gav04}.
Our model highlights two particular outcomes of hybridization which
can be of general importance. The first concerns the evolution of
Mullerian mimicry when a common hybrid coloration pattern spreads
across a system of sympatric or parapatric species as a result of
hybridization. In our simulations, this outcome corresponds to the
loss of $Z_{2}$ chromosome. The second is that the deleterious gene
flow resulting from hybridization can prevent the spread of an ecomorph
into a suitable area and, thus, can limit the species range. In our
simulations, this outcome corresponds to the formation of a tension
zone separating two specialist ecotypes at a border of the central
area. The question of whether deleterious gene flow can restrict species
ranges is currently of great interest (e.g. \citealt{kir97,cas00}).
In contrast to current discussions, in our model the source of the
deleterious gene flow is not the central populations of the same species
but rather a different species or ecotype occupying adjacent areas.
Our model was intentionally tailored for a particular case study.
Some features of the Heliconius system are probably unique (e.g. specific
genetics of the color patterns). Others may not be very common in
nature (e.g. the fact that color patterns play an important role simultaneously
in survival and mating \citep{jig01,mal89,mav06}, and thus represent
{}``magic traits'' in the terminology of \citealt{gav04}). At the
same time, our model points to some general features of hybrid speciation,
i.e. speciation via introgressive hybridization with hybrid traits
directly and significantly contributing to the survival and reproductive
isolation. Hybrid speciation can be triggered by a number of factors
increasing the likelihood of hybridization and can be driven by strong
ecological selection promoting the spread of a particular advantageous
combination of hybrid traits. The model suggests that stable sympatric
coexistence of the hybrid and parental forms should not be generally
expected. Indeed, if a particular hybrid combination of traits is
advantageous, it is expected to spread after a secondary contact between
a hybrid and parental species (if hybridization between them is possible).
In other words, the same process that leads to the creation of a hybrid
species may prevent its stable sympatric coexistence with the parental
forms. The model also shows that equal representation of parental
neutral markers in the hybrid species in unlikely. In fact, even parental
traits of ecological importance can be inherited in a hybrid species
in a relatively random fashion. The model highlights the importance
of pre-existing assortative mating and habitat segregation between
hybrids and at least one of the parental forms in simplifying the
conditions for hybrid speciation. It also shows that recombinational
hybrid speciation is not an instantaneous process but rather can take
hundreds and thousands of generations.
In contrast to many other speciation scenarios \citep{gav04}, theoretical
work on hybrid speciation has been very limited (\citealt{mcc95,bue00}
and this paper) so it is difficult to identify the most important
evolutionary factors and forces controlling its dynamics. Still existing
work suggests that the questions of spatial heterogeneity in selection,
coexistence of hybrid species with their parental species, and of
the stability of species ranges are of great importance for developing
an adequate theory. Currently, mathematical models support the belief
that hybrid speciation in animals is plausible under certain conditions.
However much more work (both empirical and theoretical) is necessary
to be able to make more definite conclusions on its importance in
nature.