Stable Heliconius butterfly hybrid zones are correlated with a local rainfall peak at the edge of the Amazon basin

Multilocus clines between Müllerian mimetic races of Heliconius butterflies provide a classic example of the maintenance of hybrid zones and their importance in speciation. Concordant hybrid zones in the mimics Heliconius erato and H. melpomene in northern Peru were carefully documented in the 1980s, and this prior work now permits a historical analysis of the movement or stasis of the zones. Previous work predicted that these zones might be moving toward the Andes due to selective asymmetry. Extensive deforestation and climate change might also be expected to affect the positions and widths of the hybrid zones. We show that the positions and shapes of these hybrid zones have instead remained remarkably stable between 1985 and 2012. The stability of this interaction strongly implicates continued selection, rather than neutral mixing following secondary contact. The stability of cline widths and strong linkage disequilibria (gametic correlation coefficients Rmax = 0.35–0.56 among unlinked loci) over 25 years suggest that mimetic selection pressures on each color pattern locus have remained approximately constant (s ≈ 0.13–0.40 per locus in both species). Exceptionally high levels of precipitation at the edge of the easternmost Andes may act as a population density trough for butterflies, trapping the hybrid zones at the foot of the mountains, and preventing movement. As such, our results falsify one prediction of the Pleistocene Refugium theory: That the ranges of divergent species or subspecies should be centered on regions characterized by maxima of rainfall, with hybrid zones falling in more arid regions between them.

Hybrid zones are narrow regions of phenotypic or genetic change between the ranges of parapatric, genetically distinct forms, typically representing a region of overlapping clines at multiple loci (Barton and Hewitt 1985;Hewitt 1988). Stable hybrid zones can be maintained by exogenous selection, that is, the fitness of genotypes is dependent on the abiotic or biotic environment and varies geographically (Endler 1977). However, hybrid zones maintained by a balance between dispersal and endogenous selection against hybrids, termed "tension zones," are not tied to fixed features in the environment (Key 1968; Barton and Hewitt 1985). One of the most interesting features of tension zones is that they have an inherent tendency to move if there are any asymmetries of selection, dispersal, or density (Bazykin 1969;Barton 1979;Dasmahapatra et al. 2002;Barton and Turelli 2011). A moving hybrid zone might easily become trapped in regions of low population density (Barton 1979; Barton and Turelli 2011). Nonetheless, in a number of empirical examples, hybrid zones have shown substantial geographic movements over periods of decades (Buggs 2007) especially in birds (Secondi et al. 2003;Krosby and Rohwer 2010;Carling and Zuckerberg 2011;Taylor et al. 2014). The stability and movement of tension zones remains a poorly explored topic in spite of its potential importance in understanding speciation and biogeography. In addition, problems caused by many invasive species also involve hybridization; thus, a greater understanding of moving hybrid zones could have important implications for conservation of biodiversity (Buggs 2007).
We studied hybridizing races of Heliconius butterflies in northern Peru. Heliconius butterflies are warningly colored Neotropical butterflies that participate in Müllerian mimicry rings, where species with convergent color patterns share the cost of educating predators about their unpalatability (Brown 1981). Heliconius erato and H. melpomene are widely distributed co-mimics that display parallel geographic variation in color pattern throughout their range, with almost 30 subspecies each currently recognized (Brown 1979;Mallet 1993;Lamas 2004). Races are usually monomorphic within their range but form hybrid zones where they meet. In the Río Mayo and upper Río Huallaga valleys in northern Peru, H. erato favorinus and its co-mimic H. melpomene amaryllis exhibit a crimson band on the forewings (known as the "postman" pattern), as well as a yellow hindwing bar (Fig. 1). In the adjacent Amazonian lowlands, the "dennis-ray" patterned H. erato emma and its co-mimic H. melpomene aglaope are characterized by a yellow band on the forewing and orange-red basal forewing patch (dennis) and hindwing rays. A narrow (ß10 km) wide hybrid zone separates the Amazonian and Andean valley forms within each species. These hybrid zones were previously documented by a study in 1984-1987(Mallet et al. 1990. Polymorphisms in warning color and Müllerian mimicry are expected to cause frequency-dependent selection against the rarest phenotype: Rarer forms will be less well recognized as unpalatable and are more likely to be attacked. Therefore, hybrid zones between divergent Müllerian mimetic races of Heliconius can be maintained by frequency-dependent selection against rare forms (Mallet 1986;Mallet and Barton 1989a,b). As with other tension zones, asymmetries in selection or dispersal between mimicry races may cause such zones to move (Mallet and Barton 1989a,b;Sasaki et al. 2002;Kawaguchi and Sasaki 2006). However, clines at dominant loci under frequency-dependent selection are liable to move even in the absence of asymmetries in dispersal or selection. For example, in the center of a cline, where allele frequencies are equal (P = 0.5), the recessive phenotype, assuming random mating, will be rarer (P 2 = 0.25) and is selected against. Consequently, a dominant allele should expand at the expense of the recessive allele because of a selective asymmetry generated entirely by genetic dominance; a process termed "dominance drive" (Mallet 1986;Blum 2002). This inherent instability coupled with the strong likelihood of some selective asymmetry between different warning patterns therefore implies that hybrid zones maintained by frequency-dependent selection should almost always be mobile. Indeed, a hybrid zone between Panamanian variants of the postman pattern in H. erato has moved ß40 km in 20 years (Blum 2002).
Variants of the dennis-ray pattern are found throughout the Amazon basin. In contrast, postman pattern races are found in disjunct regions around the periphery of the Amazon, from Central America to southeastern Brazil. Mallet (1993) suggested that these curious disjunctions might be due to the rayed pattern originating in the Amazon and spreading outwards at the expense of the ancestral, postman pattern. Recently, DNA sequence data from the optix color patterning locus ) supported this hypothesis: optix alleles specifying dennis-ray patterns are monophyletic across multiple races within both H. erato and H. melpomene , and younger than the corresponding postman races, even though the same patterns arose independently in these two mimetic species (Supple et al. 2013).
It is unclear why dennis-ray phenotypes replaced the postman phenotypes. One possibility is that rayed individuals are selectively favored because they are protected by a large preexisting dennis-ray mimicry ring found predominantly in the Amazon basin (e.g., H. xanthocles, Neruda aoede among others). In comparison, postman mimicry is more or less restricted to H. erato and H. melpomene and two other species with restricted ranges (H. besckei and H. timareta -Mérot et al. 2013). Alternatively, the dennis-ray phenotype could be more effective at teaching predators that the bearer is unpalatable, or superior in lowland Amazonian conditions (Mallet 1993). Another hypothesis is dominance drive, because all but one of the alleles determining the dennisray phenotype are dominant, both within H. erato and within H. melpomene (Mallet 1989).
Therefore, dominance and mimicry advantages for rayed phenotypes predict that the hybrid zone should be moving westwards toward the Andes. According to predictions in Mallet et al. (1990, p. 932), the movement would be expected to have been around 12.5 km in the 25 years elapsed since the previous study. In addition, the study area has undergone considerable deforestation and probable climate change. Habitat alteration could affect the hybrid zone in a number of ways. If different phenotypes are favored in different habitats (e.g., closed forest vs. more open areas; Blum 2008), deforestation might cause the more open-adapted race to increase at the other's expense. Alternatively, habitat loss might reduce predator populations and consequently relax selection, allowing a broadening of the hybrid zone. Increasing temperatures as a result of global climate change might also lead to changes in hybrid zone shape and position (Taylor et al. 2014). For example, warming temperatures might effectively "lower" the mountain passes (1000 m) separating the Río Mayo and upper Río Huallaga valleys from the Amazonian lowlands and allow increased gene flow between the two.
This Peruvian system of hybrid zones has become an important hybrid zone example. Many subsequent molecular studies of Heliconius and other mimetic butterflies have been carried out in the Mayo-Huallaga region of Peru Joron et al. 2001Joron et al. , 2006Joron et al. , 2011Baxter et al. 2010;Counterman et al. 2010;Dasmahapatra et al. 2010;Hines et al. 2011;Heliconius Genome Consortium 2012;Nadeau et al. 2012;Martin et al. 2013). In view of this near model status and the span of time since the original study (ß1986), it seemed worth resampling this hybrid zone system to test for changes over the last 25 years. Specifically, we address the following questions. (1) Have the positions of the hybrid zones moved, as expected due to dominance drive, selective asymmetries, and environmental change? (2) If they have not, what factors might be responsible? We use simulations to show how environmental variation might halt a moving cline. Previously we hypothesized that the hybrid zones are associated with a region of exceptionally high rainfall that could cause a population density trough for sunshine-dependent butterflies and trap an otherwise mobile hybrid zone (Mallet et al. 1990;Mallet 1993). We test this hypothesis explicitly using modern GIS data. (3) Have the widths of the hybrid zones and the strength of natural selection changed over 25 years, as expected due to environmental change? The previous study estimated an averaged selection pressure across all of the loci, but for the first time, we here analyze linkage disequilibria across the hybrid zones to estimate absolute selection pressures on each of the multiple loci controlling color pattern variation. Mallet et al. (1990) Fig. 1) crossing the hybrid zone at approximately right angles. This transect was near to the majority of collecting sites. The positions of collecting sites along the transect were determined by dropping a perpendicular line from each site onto the transect and calculating the distance from the intersection to the transect's start. We followed Mallet et al.'s (1990) approach but recorded the position of individual specimens along the transect, before assigning specimens to 1 km interval bins along the transect. Wings were removed and used to score genotypes of color pattern loci. The remaining tissue was preserved in NaCl-saturated Dimethyl sulfoxide (DMSO) at −20°C. Wings and preserved tissue samples are stored in the Department of Biology at the University of York, UK.

GENETICS OF COLOR PATTERNS
A detailed understanding of the Mendelian genetics of color patterns in the Amazonian and Andean races of H. erato and H. melpomene (Mallet 1989;Mallet et al. 1990) allows specimens to be genotyped using their color pattern phenotypes. We briefly summarize the color pattern genetics here. See also a fuller description in the Supporting Information (Supplementary Methods: S1).
In H. erato, three unlinked loci determine the major color pattern elements. The codominant D locus determines the presence/absence of proximal forewing patch known as "dennis" as well as hindwing rays, and the yellow versus red color of the forewing band. A dominant locus Sd determines the width of the forewing band. Sd and the dominant Cr locus interact to produce the yellow hindwing bar of the postman pattern race.
In H. melpomene, four dominant loci determine color pattern differences. The D locus determines the presence or absence of the dennis and rays. The Yb locus determines the presence or absence of the yellow hindwing bar. The N locus interacts with another locus, B to determine the shape and color of the forewing band. The N and Yb loci are tightly linked (with recombination rate, c 0.005), and the B and D loci are thought to be moderately (c 0.135) or tightly linked (Sheppard et al. 1985;Mallet 1989;Baxter et al. 2010). Amazonian alleles are dominant in both species, except for the B locus which is dominant in the Andes. The H. melpomene N locus can be difficult to score due to epistatic effects of another locus, M (Mallet 1989, Supplementary Methods S1). M has never been mapped but is almost certainly linked to the B-D region (S. Baxter, pers. comm.), and its existence may explain curious features of (apparent) disequilibria measured among mapped loci (see Results).

ITERATIONS OF THEORETICAL CLINES
Theoretical allele frequency clines were generated by iterating deterministic models of frequency-dependent selection acting in three-and four-locus hybrid zones, following Mallet and Barton (1989a,b). Multilocus iterations were used, as clines for dominant loci are asymmetric and therefore cannot be modeled appropriately by symmetric sigmoidal curves. Moreover, strong selection and linkage disequilibria between dominant and codominant loci distort cline shape away from weak selection, single-locus analytical theory, even for codominant loci (Mallet and Barton 1989a,b;Mallet et al. 1990). The R code (R Core Team 2013) used is made available in the Supporting Information, along with additional details of practical implementation (Supplementary Methods: S2). Our model of selection reduced the frequency of a genotype in linear proportion to its phenotypic rarity. The model for relative fitness (W) for each locus, from equation (4) in Mallet and Barton (1989a,b) was: ( 1 a ) Here, s is the frequency-dependent selective coefficient, P Aa etc. represent the frequencies of genotypes, and h is the dominance of allele A, (h = 0.5 if codominant, or h = 1 if allele A is completely dominant to a). Per locus fitness were then combined multiplicatively across loci. For H. erato, we iterated the model with three unlinked loci: one codominant locus (D) and two loci with complete dominance (Cr and Sd). For H. melpomene, we iterated the model for four dominant loci (B, D, N, and Yb), with moderate linkage (c = 0.135) between B and D and tight linkage (c = 0.005) between N and Yb (Mallet et al. 1990).
For each species, we generated a single iterated equilibrium solution with cline shapes for each of the loci. We then used likelihood to independently fit these theoretical cline shapes to the observed data. We tested null hypotheses that cline widths and positions remained unchanged between 1986 and 2011 by comparing summed likelihoods of clines fitted to the two separate datasets with those when clines were fitted to the combined dataset. Test statistics were calculated as twice the log-likelihood ratio (−2 L), which approximates a χ 2 distribution, with degrees of freedom equal to the extra numbers of parameters added.
To test whether a cline that is moving due to dominance drive may be halted (Barton 1979;Mallet 1993;Barton and Turelli 2011), we ran iterations in which a moving cline encounters (1) a region of low population density and (2) a separate, unlinked allelic cline tied to an environmental boundary. For details on how (1) and (2) where modeled, see the Supplementary Methods (S3).

ESTIMATES OF LINKAGE DISEQUILIBRIA
We used the 2011 dataset to estimate linkage disequilibria, thereby allowing comparison with the 1986 dataset. For polymorphic collecting sites, we calculated the maximum likelihood estimate of the linkage disequilibrium parameter D assuming Hardy-Weinberg equilibrium (Hill 1974) for all pairwise comparisons of loci. The maximum likelihood for D was compared with the likelihood when D = 0 using a likelihood ratio tests. To place these measures of disequilibrium between −1 and +1, we expressed them in terms of the correlation coefficient In H. erato, all loci are unlinked, so there is no expectation that either gene frequencies or pairwise disequilibria should be heterogeneous within populations. We obtained the maximum likelihood estimate of the correlation coefficient R AB over all three genes. In theory, linkage disequilibria are expected to peak in the center of a hybrid zone and decline to zero in the tails (Barton 1982;Mallet and Barton 1989a,b), but in practice there is often considerable heterogeneity between sites (Mallet et al. 1990;Szymura and Barton 1991). To estimate the maximum linkage disequilibria between loci, we combined the 1986 and 2011 datasets, and regressed the correlation coefficient R against the product of mean allele frequencies across loci, pq, with the constraint that the intercept of the model (when pq = 0) is equal to zero. The predicted value for R when pq = 0.25 is interpreted as the maximum value for linkage disequilibria in the hybrid zone (Szymura and Barton 1991), with confidence intervals estimated from the model.

ESTIMATES OF SELECTION PRESSURE
The width, w A , of a cline in allele frequency at a single-locus A is defined as the inverse of the maximum slope. Assuming weak selection, the equilibrium width, w A , will be proportional to the ratio of migration, σ (measured as the standard deviation of parent-offspring distances along a linear axis) to selection, s A : where K A is a constant depending on the type of selection (Barton and Gale 1993). Under frequency-dependent selection, K A = Ý8 for a dominant gene and K A = Ý12 for a codominant gene (Mallet and Barton 1989a,b). In other words, frequency-dependent selection is more effective for dominant genes, leading to narrower clines for the same s A . However, for stronger selection pressures, K A is no longer constant, but increases with selection pressure (Mallet and Barton 1989a,b). Nonetheless, the ratio of K A between dominant and codominant genes is approximately constant (Ý8/Ý12), even at high selection pressures (Mallet and Barton 1989a,b: Fig. 2). This result can be used to estimate relative selection from relative cline widths. Assuming weak selection, the linkage disequilibrium, here measured as the gametic correlation coefficient, R AB , between two loci generated by dispersal in a hybrid zone is approximately proportional to the square of the dispersal rate (σ 2 ) divided by the product of recombination rate (c AB ) and widths (w) of the clines (Barton and Gale 1993, eq. (1); Mallet and Barton 1989a, eq. (9)): Given estimates of recombination and cline width, equation (3) can be solved to provide an estimate of average selection per locus, √ s A √ s B . For larger selection pressures and/or cases involving more than two loci, computer iterations can be used to estimate average selection per locus given a value of average linkage disequilibria (Mallet and Barton 1989a,b) . Relative widths of the clines then enable estimates of the relative and absolute selection pressures on each locus. For instance, in H. erato suppose the averaged selection pressure as inferred from linkage disequilibria and itertions is s. Given equation (3), one can assume that for three loci, A, B, and C, under strong selection that: and that: We would like to estimate individual per locus selection pressures s A , s B , and s C from this information. Therefore, we can estimate relative selection pressures k and j from relative differences in cline width multiplied by the ratio of constants K i from equation (3) via a simple pair of simultaneous equations.

CORRELATIONS WITH CLIMATE
Using overall estimates of cline shape from combined 1986 and 2011 datasets, we tested the hypothesis that the position of the hybrid zone in H. erato is correlated with rainfall, as hypothesized by Mallet et al. (1990). Many climatic variables could have been used, but the amount of rainfall, which is inversely correlated with insolation and temperature in this tropical climate, is expected to be a major effect on the levels of activity of butterflies. Heliconius erato was used rather than H. melpomene because more samples were available, and the D locus was used as heterozygotes at this locus can be reliably identified. However, because all clines have similar centers (Tables 1 and 2), the results hold for all loci in both species. We obtained 30 arcsec (ß1 km resolution) gridded data for annual precipitation for the study area from the Worldclim website (http://www.worldclim.org; [Hijmans et al. 2005]). We pooled the butterfly data into the cells of the grid, with the cells then comprising sample units. We used a generalized linear model (GLM) with binomial errors and logit link to model the relationship between the proportion of heterozygotes and annual precipitation. We detected overdispersion in the model and so corrected the standard errors using a quasi-GLM model where the variance is given by φμ, where μ is the mean and φ is the dispersion parameter. To validate the model, we plotted the Pearson's residuals against fitted values and the explanatory variables, and tested them for spatial autocorrelation, which can violate the assumption of independence of errors, using a Mantel test applied to Pearson's correlation coefficient. In addition to this gridded  (Table 1). For dominant clines, small samples on the Amazonian side of the clines often failed to detect any recessive homozygotes, which are necessary for a nonzero estimate of recessive allele frequency following Hardy-Weinberg principles. Thus, apparent deviations from the model result from small sample sizes, not mis-specification of the model. climate dataset, we also directly obtained data on rainfall for towns in the study area from the website http://www.agteca.com. Table 1 compares the likelihoods and parameters for the bestfitting theoretical clines estimated for the combined and separate time periods. These are plotted with the data in Figure 2A, B. Fitting separate clines to data from 1986 and 2011 resulted in a significant increase in likelihood over the combined datasets only for the H. melpomene B locus. All other loci were adequately described by single clines fitted to the combined datasets. To test whether the increase in likelihood for the H. melpomene B locus was due to changes in cline position or slope, we fitted separate clines to each time period but allowed only centers or only slopes to vary. Specifying different slopes with a single center gave a significant increase in likelihood (the change in ln likelihood, L = 3.70, 1 df, P = 0.007) compared to the combined dataset. Specifying different centers with a single slope gave no significant increase ( L = 0.40, n.s.).

SPATIAL AND TEMPORAL CHANGE
These results suggest that over approximately 100 generations since 1986, the clines remained effectively stationary to an accuracy of ±1 km (this is approximately 1/3 of the per generation dispersal distance in H. erato, or 1/6 that in H. melpomene-see estimates below). Among possible reasons for the lack of cline movement, we found that either a reduction in density or interaction with an unlinked, exogenously selected cline tied to the environment can both halt clines at loci that are moving due to dominance drive (Fig. 3). Which of these two alternatives is more powerful depends on the gradient and length of the slope in increasing density faced by the moving cline as it exits the density trough, versus the strength of selection on the unlinked cline.
Clines that are the result of secondary contact between neutral traits are expected to decay (Endler 1977). The width of a neutral cline t generations after the populations met in an abrupt step

Likelihoods and parameter estimation for H. erato color pattern locus clines. Support limits in brackets. Cline centers and widths are not significantly different between
1986 and 2011 ( Fig. 2A) Likelihoods and parameter estimation for H. melpomene color pattern locus clines. Support limits in brackets. Cline centers and widths are not significantly different between 1986 and 2011 except for the B locus width (Fig. 2B).

Log-
where σ is the parent-offspring dispersal distance (Barton and Gale 1993). Assuming σ 2-3 km (Mallet et al. 1990, and see below), with four generations per year and an initial width of 7.6-16.7 km in 1986 (depending on the locus), if the color pattern clines were the product of neutral mixing following secondary contact they would be expected to have widened to between 57.8 and 92.0 km by 2011. The widening would have been readily detected, as it is well outside the support limits for the widths measured in 1990 (Tables 1 and 2). From the lack of significant broadening of all clines in both species, we can infer the color patterns are under selection, and that similar selection pressures must have maintained the hybrid zone in 1990 and 2011. No cline shows evidence of movement, and only locus B in H. melpomene shows a significant change in cline width. However, rather than broadening it has narrowed somewhat, suggesting, if anything, increased selection.

ESTIMATES OF OVERALL SELECTION IN H. ERATO
In the 2011 H. erato data, significant (P < 0.05) linkage disequilibria were found only between D and Cr loci, at two sites on the Andean side of the hybrid zone (Table 3). Joint disequilibria were estimated at only at four sites polymorphic for all three loci. As in the better sampled data from 1986, disequilibria were mostly positive, but only one site was significantly different from 0 (Table 3), with R = 0.22. When combined with disequilibria estimated from all three loci in 1986 (Mallet et al. 1990), regression analysis gave an estimated maximum R = 0.35 (support limits: 0.20-0.49; Supporting Information Fig. S1). This suggests an average per locus selection s = 0.22 (0.13-0.31; see Fig. 2B of Mallet and Barton 1989a,b), hardly different from that inferred for the 1986 data alone (s = 0.23, Mallet et al. 1990).

SELECTION PRESSURES ON EACH MIMICRY LOCUS IN H. ERATO
The D cline was 0.81 (0.72, 0.92) times the width of the Cr cline, and 0.77 (0.68, 0.87) times the width of the Sd cline (support limits were generated by increasing/decreasing the widths of clines until the log likelihood of the fits to the combined data dropped by 2 units). Given that Cr and Sd are dominant whereas D is codominant, selection pressures on Cr and Sd relative to D are k = 0.44 (0.35, 0.56) and j = 0.40 (0.31, 0.50), respectively (from eqs. (4), (5)). Estimates of absolute selection pressures on each locus are therefore s D 0.38, s Cr 0.17, and s Sd 0.15. (Note: these estimates are based on interactions between loci, so levels of uncertainty are given only for the average selection, s, and relative selection pressures, k, j. It is not easy to infer confidence intervals for the absolute selection pressures on each locus, but the support intervals for k, j give a good idea of the uncertainty. These relative selection pressures are as might be expected under selection via visual predation: The D locus affects far more of the wing surface than do Cr or Sd. Although we cite separate selection pressures for Cr and Sd, they are not significantly different, whereas selection on the D locus is significantly stronger than on the other two, as is readily seen from the significant differences in cline width (Table 1). These separate per locus selection values are similar to those estimated based on a mark-recapture study (s D 0.33, s Cr 0.15, and s Sd 0.15, ). Dispersal distance can be calculated using the average cline width of 9.24 km. The average selection pressure of s = 0.22 implies cline width of about 3.8-fold of dispersal (Mallet and Barton 1989a,b: Fig. 2), so dispersal, σ 2.6 km per generation.

DISPERSAL IN H. MELPOMENE
Disequilibria were higher and more often significant in H. melpomene; some sites had significant disequilibria for every pairwise comparison of loci (Table 4). As expected, pairwise disequilibria are higher for linked loci (B vs. D, N vs. Yb), but even unlinked loci showed much stronger disequilibria than observed in H. erato, in spite of smaller sample sizes. The recombinant genotype N-ybyb was over twice as common in the hybrid zone as the reciprocal recombinant, nn Yb-(24:10 in combined data, significantly different from the 1:1 expected, L = 2.97, 1 df, P = 0.015). This imbalance is almost certainly due to interaction with the unscored M locus (Supplementary Methods: S1): Some nn ybyb mm genotypes were probably mis-scored as recombinant N-ybyb (Mallet et al. 1990). Therefore, in H. melpomene we base our disequilibrium estimates between unlinked genes on disequilibria only between B or D and Yb, which is not prone to hidden epistasis.
Regression analysis gave an estimate of maximum R = 0.56 (0.42, 0.71) for unlinked loci D and Yb (Fig. S2). Reading from Figure 2B in Mallet and Barton (1989a,b), this suggests average per locus selection of 0.31 (0.22, 0.43). Linkage disequilibria for the also-unlinked loci B and Yb while not independent of the D-Yb disequilibrium estimates, due to physical linkage between D and B, are expected to be lower due to the greater width of the B cline (see eq. (3)). Likelihood fits show that widths of D, N, and Yb clines are not significantly different, but that all of them are 0.70 (0.63, 0.78) times as wide as the B cline. The width difference is presumably caused by weaker selection on B, and the b-Yb disequilibria are indeed smaller (although not significantly so). Regression analysis gives the maximum disequilibrium R = 0.45 (0.22, 0.68) (Fig. S2), equivalent to a combined selection pressure of √ s B √ s Yb = 0.24 (0.11, 0.40). Given that the cline widths of N, Yb, and D do not differ significantly, then s Yb 0.31 (from the disequilibria of D and Yb above), suggesting a value  For the linked loci B-D, as well as N-Yb, the simple regression analysis used here fails to give reasonable values of maximum disequilibrium (R 1) because it does not accommodate the possibility that R can be no higher than one (Figs. S3 and S4). However, these gametic correlations are indeed close to 100% across the hybrid zone (Table 4), as expected for tightly linked loci under selection across a narrow hybrid zone. Probably, disequilibria between N and Yb are even higher than shown in Figure S3 because scoring difficulties for N lead to overestimates of numbers of N-ybyb recombinant genotypes (see M locus discussion above). With R 1 for N and Yb, the two loci become so completely intertwined that they act largely as a single locus

Maximum likelihood estimates of linkage disequilibrium, R, for H. erato, 2011. Significant values shown in bold.
Allele frequencies 19 and s D ≈ 0.31) is somewhat compromised due to disequilibria: Differences in selection between B and D are therefore liable to be underestimated, whereas the per locus selection based on the assumption of lack of linkage on each locus is liable to be overestimated because each locus compounds the selection on the other via the strong gametic correlations. Bearing these difficulties in mind, however, our estimates of selection were again similar to earlier averaged estimates for H. melpomene (s Yb = s N = s D = s B 0.25, Mallet et al. 1990).
Dispersal in H. melpomene is most easily inferred via cline widths and disequilibria between unlinked loci Yb and D (N is ignored due to scoring difficulties). These have conveniently similar cline widths (Table 2), and therefore similar selection pressures (s 0.31). Given three dominant loci, Figure 2 of Mallet and Barton (1989a,b) shows that s = 0.31 corresponds to cline width 1.89 times the dispersal distance. The cline widths of 10.9 km on average therefore imply σ 5.8 km.
Both selection and gene flow are here estimated in H. melpomene to be rather greater than earlier (Mallet et al. 1990), but this is mainly due to differences in methods of estimating maximum R, which led to somewhat lower previous estimates.
The actual values of disequilibria are similar in the two time periods. Heliconius melpomene probably experiences somewhat more selection as a result of its higher dispersal across the hybrid zone, so that a higher proportion of H. melpomene than H. erato is exposed to predators trained to warning colors different from those in their parental populations.

CORRELATIONS WITH RAINFALL
We found a strongly positive correlation (P < 0.001) between the proportion of heterozygotes at the H. erato D locus and annual precipitation (in meters) with an intercept of zero (−8.04 ± 0.96 on a logistic scale, or 0.00 as a proportion) and a positive slope (2.91 ± 0.43 on a logistic scale, or 0.95 as a proportion). The correlation explained 42% of the deviance. The model fit is shown in Figure S5. Figure S6 plots monthly rainfall across the year for five key sites along the transect, showing that even the driest month at the center of the hybrid zone is wetter than the wettest month in Tarapoto on the Andean side of the hybrid zone. These data also show that the rainfall peak in the center of the hybrid zone transect at Pongo del Caynarachi, averaged across an 18-year period, was underestimated by the interpolated gridded climate data used in the correlation analysis. A map of interpolated rainfall is shown in Figure S7. Otherwise, the values from the gridded data were similar to the weather station data. A Mantel test applied to the scaled Pearson residuals revealed no significant spatial autocorrelation (10,000 permutations, r = −0.07, P = 0.85). Along the transect, rainfall peaked at the base of the Andes, corresponding to Pongo del Caynarachi at the center of the hybrid zone system (Fig. 3).

Discussion
Previous theoretical and empirical work has demonstrated that Heliconius color pattern clines can and do move across the landscape (Mallet 1986;Mallet and Barton 1989a,b;Blum 2002;Sasaki et al. 2002;Kawaguchi and Sasaki 2006). Here, hybrid zones that had been predicted to move were shown to have remained remarkably stationary in a zone of exceptionally high rainfall over a period of 25 years. Mallet et al. (1990) and Mallet (1993) hypothesized that predominantly cloudy conditions and heavy rainfall at the base of the Andes could reduce recruitment and increase mortality, and hence lower equilibrium butterfly densities. Regions of low density can prevent cline movement even when the advancing race has a significant fitness advantage (Barton 1979;Sasaki et al. 2002;Barton and Turelli 2011), and demographic stochasticity does not strongly affect the outcome (Kawaguchi and Sasaki 2006;Barton and Turelli 2011). Figure 3 shows a theoretical example where a density trough halts movement of a frequency-dependent cline advancing due to dominance drive. The advancing dominant cline is halted because, as it attempts to climb out of the density trough, an increasing density gradient causes asymmetric migration in favor of the recessive allele. This can prevent the dominant phenotype from reaching the critical phenotypic frequency of 50% in populations the other side of the trough, as is required to allow the clines movement to continue. Figure 3 also shows that, just before being halted, the cline briefly accelerates into the region of low density because at this stage the decreasing density gradient and dominance together promote movement in the same direction. A moving cline can also become trapped when it encounters and interacts with stationary clines maintained by exogenous selection along an ecotone (Fig. 3). We have previously shown that different tension zones, some of which may be mobile and some stationary, will also be attracted together, and that after contact they will move along together (Mallet and Barton 1989a,b;Mallet and Turner 1998: Fig. 16:3). As far as we know, the potential attraction of tension zone clines to other clines has not previously been discussed, but may be important in understanding geographic localization of hybrid zones in general. Cline attraction is here due to linkage disequilibrium: Moving clines become subject to increased per locus selection as they cross the ecotone due to disequilibrium with clines under exogenous selection. Provided selection and thus disequilibria are strong enough, this can counter the selective asymmetry that causes movement. The attraction between clines again causes a slight acceleration just before the moving cline is halted (Fig. 3).
Stability of a mimicry cline might occur if the second locus coded for adaptation to a host plant or local climatic conditions associated with altitudinal change or aridity. Heliconius himera may provide an example. This species hybridizes with H. erato at low frequency but forms narrow hybrid zones wherever the two species meet; H. himera appears limited to gallery forest in arid inter-Andean valleys (Fig. 1;Mallet 1993;Jiggins et al. 1996). Alternatively, color pattern clines themselves may also be under exogenous, nonmimicry associated selection, with the rayed form a direct adaptation to Amazonian rainforest (Mallet 2010).
We tentatively suggest that ecotone-based scenarios seem less likely than the rainfall density trough hypothesis in the Peruvian system studied here. First, there is little evidence for genomic differentiation between the races other than at color pattern loci Nadeau et al. 2012;Martin et al. 2013). Second, recent collections from the Alto Río Mayo watershed in the northern part of the postman phenotype's range have revealed the same H. melpomene and H. erato hybrid zones in the Andes at ß1200 m near the headwaters of the Río Mayo, approximately 60 km from the nearest Amazonian lowlands ( Fig. 1; M. Joron, pers. comm.). Pure rayed H. erato emma and H. melpomene aglaope also occur at elevations up to 1500 m in the Río Pozuzo valley some 400 km to the south of Tarapoto (Mallet 1993). These distributional patterns therefore suggest that the dennis-ray phenotype is not itself an adaptation to lowland habitats or climate. These geographic patterns also seem inconsistent with trapping by interactions with other clines adapted to local conditions. Indeed, the upper Río Mayo location of the same hybrid zones is particularly interesting because it is not so obviously associated with ecotones or climatic features, raising the intriguing possibility that it may be free to move. Some previous studies of H. erato have concluded that hybrid zones associate with ecological adaptation leading to parapatric speciation to distinct habitats (Jiggins et al. 1996;Arias et al. 2008;Blum 2008). However, wing pattern differences do not necessarily constitute adaptations to the local environment. Instead, mobile tension zones may simply become trapped by density troughs or clines at other, ecology-adapted loci. A potential example of this may be found along the Amazon river, which marks a region of major Heliconius subspecies turnover, despite the fact that butterflies do disperse across it (Rosser et al. 2012). Once clines have stabilized, there is then potential for local adaptation and genomic divergence at linked loci, ultimately leading to speciation; the socalled "islands of divergence" model of speciation (Feder et al. 2012).
We found that the widths of the clines, as well as their positions and linkage disequilibria among loci were remarkably similar between 1986 and 2011. This provides evidence that relatively constant selection against foreign phenotypes maintains these color pattern clines. Significant deforestation over the last 25 years on the Amazonian side of the cline might have been expected to have reduced predator populations and consequently selection, giving rise to a potential widening of the cline. Jacamars and tyrant flycatchers are thought to be the most important visual predators responsible for the frequency-dependent selection on Heliconius (Chai 1986(Chai , 1996Pinheiro 1997;Langham 2004;Merrill et al. 2012). However, like Heliconius, jacamars and flycatchers are often gap specialists and may therefore respond in similar ways to deforestation, thus approximately maintaining the relative frequencies of predators and prey to yield similar levels of selection.
There has often been an implicit assumption that hybrid zones are merely stationary markers of secondary contact between formerly isolated populations (Mayr 1963;Haffer 1969). Tension zone theory shows that while zones of hybridization might usually stabilize in shape as a result of a balance between selection and gene flow, they are highly likely to be positionally unstable, including in specific models of mimicry hybrid zones (Barton and Hewitt 1989;Mallet and Barton 1989a,b;Sasaki et al. 2002;Kawaguchi and Sasaki 2006;Barton and Turelli 2011). A large number of moving hybrid zones are now known (Buggs 2007), with many good examples in vertebrates (Arntzen and Wallis 1991;Bronson et al. 2003;Secondi et al. 2003;Haas and Brodin 2005;Carling and Zuckerberg 2011) and a few from invertebrates (Britch et al. 2001;Dasmahapatra et al. 2010). Other hybrid zones may be positioned as a result of stationary ecological adaptations to physical ecotones, but even here, global environmental shifts such as climate change is liable to alter the null point of adaptation between two differently adaptive races; in one case, a bird hybrid zone has been suggested to have moved in response to climate change (Taylor et al. 2014).
Thus, the assumption that hybrid zones are stationary relics of former secondary contact must be abandoned. Whether secondary contact occurred, many tension zones should manifest an inherent tendency toward motion, and we fully expected this to be the case in Heliconius, especially given documented zone movement in Panama (Blum 2002). In contrast, the current article shows that two hybrid zones are in almost exactly the same positions and under the same kind of selection today as they were 25 years ago, at the Easternmost edge of the Cordillera Escalera of the Andes (Fig. 1). This is all the more remarkable as, although the hybrid zones are only about 10 km wide, they are very long. Limited samples show that the same hybrid zones are associated with the Escalera ridge from beyond Rioja, San Martín to the North, and South to the Boquerón del Padre Abad, between Tingo Maria and Pucallpa in Huánuco, a length of over 400 km. Previous reviews give the impression that most other hybrid zones are also static (Barton and Hewitt 1989), although no systematic survey of the proportion of hybrid zones showing stasis has been carried out, to our knowledge. If most are indeed static then the question remains why this should be the case.
Using a largely theoretical argument, Barton and Hewitt (1989) came to the conclusion that spatial inhomogeneities in population structure were largely responsible for trapping hybrid zones near where they were first formed. They cited the case of a grasshopper (Podisma) hybrid zone that approximately followed a density trough near the ridge top of the Pyrenees, as predicted in a two-dimensional model. The stasis of woodland bird hybrid zones in the largely treeless Great Plains of North America may have a similar explanation (Moore and Buchanan 1985;Rising 1996;Mettler and Spellman 2009). In the current study of butterflies, as well as in a number of other contact zones in Müllerian mimetic races of butterfly (Dasmahapatra et al. 2010), hybrid zones are restricted to a region characterized by catching the Westernmost orogenic rainfall of the Andes, and there is a significant correlation of the current hybrid zones with rainfall. Although we cannot entirely rule out adaptations to local conditions, it seems likely that the stasis we have observed in the hybrid zones of Heliconius and other butterflies is due to trapping by this climatic feature because rainy weather tends to reduce butterfly productivity and population density.
A longstanding hypothesis for the origin of new races and species in tropical America is the refugium theory, in which a reduction of forest cover during dry periods in the earth's history induced allopatric speciation in fragmented forest "refugia" where precipitation remained high (Haffer 1969(Haffer , 2008Brown 1979). The local correlations between high precipitation and the presence of hybrid zones we document here thus contradict one prediction of the refugium theory; that hybrid zones are expected to occur in drier or marginal areas between currently wet and putatively more stable refuge areas (Brown 1979(Brown , 1981. Conceivably, of course, these races did indeed form in allopatric refugia, but then hybrid zones between them moved to their current location following secondary contact. If so, our results suggest that current distributions may not be very useful guides to the locations of past refugia, a purpose for which Heliconius subspecies distributions have often been used (Brown 1979). Our findings therefore add to a growing body of work casting doubt on refugium theory (Endler 1977;Beven et al. 1984;Nelson et al. 1990;Bush 1994;Moritz et al. 2000;Wilf et al. 2003;Whinnett et al. 2005;Dasmahapatra et al. 2010).

Supporting Information
Additional Supporting Information may be found in the online version of this article at the publisher's website: Figure S1. Least squares regression of linkage disequilibrium (R) for H. erato estimated from all three loci against the product of average allele frequencies, using only sample sizes greater than 20. Figure S2. Least squares regression of linkage disequilibrium (R) for H. melpomene unlinked loci using sample sizes greater than 17. Figure S3. Linkage disequilibrium (R) for H. melpomene linked loci N and Yb plotted against the product of the average allele frequencies. Figure S4. Linkage disequilibrium (R) for H. melpomene linked loci b and D plotted against the product of the average allele frequencies. Figure S5. Rainfall plotted against proportions of heterozygotes at the H. erato D locus, with model predictions. Figure S6. Monthly rainfall data at sites across the hybrid zone. Figure S7. Gridded annual rainfall data overlaid onto a SRTM 90 m resolution digital elevation model. Table S1. Genotypic data for Heliconius erato. Table S2. Genotypic data for Heliconius melpomene.