Introduction

Natural hybridization between pairs of species can represent secondary contact following a past speciation event, or primary divergence during the process of speciation itself. The study of hybrid formation in zones of sympatry between two species can therefore provide insights into the evolution of reproductive isolation and the likely outcomes of its breakdown (Harrison, 1990). Hybridization rates between taxa vary widely (Ellstrand et al., 1996; Jiggins and Mallett, 2000). When hybridization is sufficiently frequent, zones of sympatry between the parental species can contain an abundance of individuals with intermediate characters, giving rise to a unimodal character distribution (for example, Burgess et al., 2005). When hybridization is instead low to moderate, such zones typically have a bimodal distribution, with few individuals with intermediate characters and many parental types (Harrison and Bogdanowicz, 1997). These characters can include morphological, physiological, genetic, or any other characters that distinguish the parental species. Note that here we use the term zone of sympatry, rather than hybrid zone, to avoid implying that hybrid individuals must be present.

It is possible for the hybridization rate to vary, not only between taxa, but also across multiple zones of sympatry of a single pair of species, as has been suggested for certain pines (Watano et al., 2004) and oaks (Williams et al., 2001). Geographic variation in hybridization rates among zones of sympatry of a single species pair could provide an opportunity to examine the effects of local conditions on reproductive isolation. If the frequency of hybrids between two species varies between zones, it is possible that reproductive isolation between them depends at least partly on local ecological circumstances, although differences in stage of speciation or of secondary contact may in principle also play a role. Identifying which ecological conditions are important to intersite differences in reproductive isolation could help us determine the relative impact of different extrinsic reproductive isolating mechanisms on reducing gene flow. For example, if two zones of sympatry that differ in hybridization rates also differ in environmental selection against hybrids, then we might conclude that postzygotic isolating mechanisms are important in isolating those species. This sort of comparison is more valuable within a species pair than between pairs of taxa with different hybridization rates because, in the latter case, between-pair differences in intrinsic isolating mechanisms confound the comparison (Coyne and Orr, 2004). For the same reasons, comparisons within a species pair might better reveal complex isolating mechanisms, such as environment-dependent behavioral isolation (Endler, 1992). Therefore, it should be of interest to quantify the extent of intersite variation in hybridization within species pairs. Despite the relatively high incidence of angiosperm species that hybridize (Ellstrand et al., 1996), such an approach has rarely been taken with flowering plants (in Gaillardia by Heywood, 1986; in Helianthus by Buerkle and Rieseberg, 2001).

We can estimate the rate of hybridization in a zone of sympatry by identifying individuals intermediate for traits that distinguish the parent species. However, the choice of traits can affect our conclusions concerning the hybridization rate. If natural selection alters the patterns of introgression of some traits, the apparent rate of hybridization in a zone of sympatry will depend on whether it was measured using neutral or non-neutral traits. Neutral markers can indicate patterns of genetic structure different from those of traits under natural selection (Martinsen et al., 2001; McKay and Latta, 2002, Streisfeld and Kohn, 2005). Theory predicts that introgression at a neutral locus will depend on its map distance from loci under selection, so clines of neutral markers that are in linkage disequilibrium with selected loci will resemble those of the selected loci, whereas unlinked neutral loci might have wider clines (Barton and Hewitt, 1985; Rieseberg et al., 2003; Gonzalez-Rodriguez et al., 2004). As a result, different neutral markers might also produce conflicting indications of the rate of hybridization.

In this study, we compared measurements of hybrid frequency based on floral morphology to the same measurements based on neutral genetic markers in two zones of sympatry between Ipomopsis aggregata and I. tenuituba (Polemoniaceae). These two closely related species (Wolf et al., 1993) were considered by Grant and Grant (1965) as a prime example of reproductive isolation because of the behavior of the hummingbird versus hawkmoth pollinators. However, the two species hybridize frequently in sympatry, suggesting some breakdown of reproductive isolation (Grant and Wilken, 1988). A monocarpic life history, showy flowers and absence of seed bank have made Ipomopsis a frequent model for study of evolutionary processes, including selection in hybrid zones and speciation (Campbell, 2004; Waser and Campbell, 2004). Studies of one Ipomopsis zone of sympatry have documented strong selection on floral traits, including color (Campbell et al., 1997; Melendez-Ackerman and Campbell, 1998), and environmentally dependent selection during the vegetative stage (Campbell and Waser, 2007), but comparatively little is known about the molecular genetic structure of any of the zones (Wolf et al., 1997; Wu and Campbell, 2005). Variation between sites in frequencies of hybrid individuals has been reported earlier (Aldridge, 2005), but that comparison was made using hybrid index scores derived only from floral morphological traits that may experience selection by pollinators (Campbell et al., 1997). It is of value to compare the genetic structure for such traits with that of neutral genetic markers.

Here, we use neutral molecular markers (RAPD) to ask if there is a variation among two zones of sympatry in frequency of Ipomopsis hybrids. We also ask if morphological and molecular markers produce similar indications of whether the distribution of phenotypes in the zone is unimodal or bimodal (hereafter, ‘modality’).

Materials and methods

Study sites and sampling

Ipomopsis aggregata (Pursh) V Grant and I. tenuituba (Rydb.) V Grant are closely related species of perennial, monocarpic herbs common in the mountains of the western United States (Grant and Wilken, 1988). Though they hybridize readily in sympatry, the species exhibit contrasting suites of floral traits consistent with the hummingbird (I. aggregata) and hawkmoth (I. tenuituba) pollination syndromes (Grant and Grant, 1965; Campbell et al., 1997) even in sympatry. Hybrids in the F1 generation are intermediate for corolla length and width and are most readily identified by a distinct orange-red corolla color (Campbell and Waser, 2007).

Two Ipomopsis zones of sympatry in western Colorado, Grizzly Ridge (GR) and Poverty Gulch (PG), display a marked contrast in the frequency of hybrids based on floral characters (Aldridge, 2005). Grizzly Ridge lies within Black Canyon of the Gunnison National Park, Montrose County, Colorado. The ridge ranges from 2375 m to 2438 m elevation, and is characterized by oak scrub and Sagebrush vegetation. Ipomopsis populations there form a mosaic of discrete, single-species patches and morphologically intermediate plants are rare despite separation of these patches by only a few meters in some cases (Figure 1a; Aldridge, 2005). Distributions of floral morphology within patches at GR closely match the distribution of morphology in reference populations of the two parent species that lie outside the zone, with only rare individuals outside of the observed range of phenotypic variation for the parent species (histograms in Figure 1a; see Aldridge, 2005 for a detailed description of the histograms). In particular, we have never observed any plants with the characteristic orange-red corolla color of F1 hybrids. Vegetation and soil analyses (unpublished data) have revealed no differences between I. aggregata and I. tenuituba patches at GR.

Figure 1
figure 1

Two Ipomopsis aggregata — I. tenuituba zones of sympatry in the Rocky Mountains of western Colorado, USA that differ in hybridization rate and spatial structure. (a) Grizzly Ridge, Montrose Co (GR); (b) Poverty Gulch, Gunnison Co (PG). Letters indicate patches of parents and hybrids sampled for floral morphology and tissue for DNA extraction. Bar plots are frequency histograms for values of a canonical discriminant function (CDF) derived using floral morphology from reference populations shown on maps (described in Aldridge, 2005). Plots correspond to groups of patches comprising either mostly parental individuals of one species or the other (both groups in panel a; A–C and K,L in panel b), mostly intermediate individuals (D–F and G–J in panel b), or to single reference populations (indicated in both panels). High values of CDF correspond to I. tenuituba morphology. Large bar plots under each site name are histograms for the entire zone; open bars indicate CDF scores for reference populations. Sample sizes are as follows: all reference populations, N=20 each; GR, N=177; PG, N=118. The arrow labeled ‘F1’ in the large bar plot for panel b indicates the mean CDF score for 33 known F1 hybrids generated by hand pollinations of parents from PG. Contour interval=20 m in both panels.

Full size image

Poverty Gulch is a tributary of the Slate River, north of Crested Butte, Gunnison County, CO, USA, approximately 110 km east of Grizzly Ridge. Elevation there ranges between 2900 m and 3231 m, and the site is predominantly montane meadow vegetation. PG is the site of an extensively studied hybrid zone (review in Campbell, 2004). The species at PG occur at opposite ends of an altitudinal cline, and are separated by a series of discrete patches with intermediate morphology, changing gradually from one species to the other (Figure 1b). These populations were first described as a hybrid zone by Grant and Wilken (1988) (see also; Campbell et al., 1997). Molecular analysis of Ipomopsis at Poverty Gulch based on one nuclear and five cytoplasmic species-specific AFLP and RFLP markers has since confirmed that all patches designated as hybrid (D-J in Figure 1b) contain some hybrid individuals with the cytotype of one species and the nuclear type of the other (Wu and Campbell, 2005). Floral morphology within the PG zone of sympatry is skewed in the direction resembling I. tenuituba. However, in marked contrast to the situation at GR, many plants here exhibit intermediate floral morphology that is outside of the range of both parental populations (Figure 1b; Campbell et al., 1997), and instead typical of the morphology of experimentally generated hybrids (Campbell et al., 2002). This greater frequency of intermediate morphologies occurs despite the greater physical distance between the two species at this site, and even between some nearest neighbor patches, compared with the situation at GR. The patches shown in Figure 1a and b represent all known patches of Ipomopsis within the areas of interest of this study, and individuals were sampled from all of the patches shown.

Morphometric data and tissue for DNA analysis were collected at both sites every year from 1999 to 2001 from flowering plants only (see Aldridge, 2005 for details on morphometric data). Because plants flower once and die, this sampling scheme prevented any resampling of individuals. However, morphometrics and tissue were not always collected from the same individuals, making direct comparison of floral morphology and genetics impossible at the level of individual plants. Plants were sampled along a transect run through the longest axis of each patch. We measured five floral traits: corolla length and width, style length and the lengths of the longest and shortest filaments. These floral morphological data were analyzed earlier (Aldridge, 2005) using canonical discriminant function (CDF) analysis; however, this method requires prior identification of morphology presumed to be typical of the two species. The CDF then calculates a linear function that maximizes the distance between the centroids of the two presumed parental types in multivariate morphological space. Our goals in the present analysis were first to compare neutral genetics to floral morphology without any a priori assumption of structure in the data, and then to include genetic information from reference populations of each species lying nearby (<2 km) but outside of each contact site. This latter procedure was analogous to using reference populations of presumed parental types to define the CDF function in analysis of floral morphology.

All morphological traits were averaged for each plant, and analyses performed on standardized (z-score) averages. Unopened buds were collected onto ice, freeze-dried within 24 h and stored at −20 °C until used for DNA extraction. Data were collected from 12 patches at each site (Figure 1a and b), and from 1 reference population for each species at each site. Sample sizes from patches within the zones (not reference populations) ranged from 10 to 30, depending on patch size.

Genetic data

DNA was extracted from bud tissue using Qiagen DNeasy plant miniprep kits (Qiagen, Valencia, CA, USA). A second extraction was performed on 10 individuals, to test for repeatability of RAPD amplification products between accessions. For inheritance tests, embryos were dissected out of seeds and extracted while fresh.

Screening for genetic markers began in 2000, using 10-mer RAPD primers (Operon Technologies, Alameda, CA, USA; kits A–N plus AB07 and AB10). RAPD was chosen for low cost of screening large numbers of primers, and ease of development. Also, AFLP screening has produced very few species-specific markers for this system (Wu and Campbell, 2005). After initial screening, 17 primers were chosen because they clearly and repeatedly produced polymorphic banding patterns in two different thermal cyclers (Eppendorff Mastercycler, Eppendorf AG, Hamburg, Germany; MJ Research PTC-100, MJ Research Inc., Reno, NV, USA), and were amplified in 425 individuals (GR=210, PG=215), producing a 48 locus RAPD phenotype for each plant.

The use of RAPD data has come under criticism for low specificity, relative ‘anonymity’ of marker bands and the possibility that multiple comigrating fragments might produce the false appearance of homologous markers (Bachmann, 1994). We minimized these problems by using only highly repeatable bands and by testing the inheritance of RAPD markers. Conspecific crosses of plants of both species were performed in 2003 on potted plants in a screenhouse (WeatherPort Co, Olathe, CO, USA) at the Rocky Mountain Biological Laboratory. The expected frequency of a dominant nuclear marker in offspring of an A_ × aa cross is either 1 or 0.5 regardless of seed parent, and the observed frequencies among 40 offsprings of each cross were compared to expectations using χ2 tests (SPSS Inc., Chicago, IL, USA).

Data analysis

We used similar methods to explore both morphological and genetic structures at the two sites. Our first goal was to test for structure without a priori assumptions of its existence, so we used maximum likelihood clustering methods to test models against our floral morphology data and RAPD phenotypes separately at each site. Both methods considered all plants together without any regard to prior identification of either the species or location within the site.

For morphology, we used the clustering method of Fraley and Raftery (2002 procedure Mclust in R; http://www.r-project.org/). This procedure tests 12 models that differ in assumptions regarding shape and volume of multidimensional clusters, for each model between 1 and 9 clusters and chooses the best-fit model based on Bayesian Information Criterion. This procedure also calculates the posterior log likelihood of the model given the data, and classifies individual observations into the clusters under the best-fit model.

To quantify the genetic structure at each site, we used the program STRUCTURE 2.2 (Falush et al., 2007). This method makes no a priori assumptions about the ancestry of test subjects, and instead clusters individuals based on their multilocus genotypes, using a Markov Chain Monte Carlo algorithm. The program calculates the probability of an individual inheriting its genotype from each of K clusters that are defined by allele frequencies, and so K can be interpreted as the number of separate parental populations represented in a sample. The version we used (2.2) explicitly incorporates the genetic ambiguity in dominant markers such as RAPDs. For each site, we used the program in two ways. First, we analyzed the 210 test subjects from GR and the 215 test subjects from PG separately to find the number of clusters (value of K) that best fit that site, by finding the posterior probabilities of K. Analyses were performed for K=1 to K=5 clusters using a burn-in of 35 000 iterations and a run-length of 35 000 iterations of the Markov chain, after test runs showed that run-length sufficient to stabilize the Markov chain and to produce highly replicable results across runs with the same parameter values. We used an admixture model with no a priori information regarding spatial population structure, to allow an individual a non-zero probability of belonging to >1 cluster.

Second, we reran STRUCTURE 2.2 incorporating prior information from the I. aggregata and I. tenuituba reference populations that were sampled outside, but near to, each of the contact sites. These reference populations (N=98 total at GR and N=97 at PG) were used as learning samples (as defined by; Pritchard et al., 2000) to define two parental populations for each site and then infer ancestry of the test individuals. Changing the degree of misclassification allowed in the learning samples (see Pritchard et al., 2000) had little effect on the outcome of STRUCTURE, and final results used v=0.05. For the test individuals, we then obtained the probability (Q) of each individual belonging to each cluster. In the case of K=2, parental individuals should have Q-values near 0 or 1 and F1 hybrids should have Q-values near 0.5. A bimodal distribution of Q indicates that a contact site has relatively few individuals of highly mixed ancestry. We assessed modality by testing frequency histograms of Q-values using the Dip test of Hartigan and Hartigan (1985 (procedure ‘dip’ in R). The dip statistic measures the maximum difference, over all sample points, between the observed distribution and the unimodal distribution that minimizes that difference and thus as dip increases, the observed distribution is less likely unimodal. The significance of dip statistics was determined using tables (‘qDiptab’) included in the Diptest package in R. A contingency analysis was used to compare the frequency distributions of Q for the test subjects from the two sites, placing individuals in bins at 0.1 intervals.

For many individuals, 90% confidence intervals on Q were very wide, and so we did not attempt to classify individuals as to whether they were parental or hybrids, but instead further evaluated averages for patches. We assessed whether the spatial genetic pattern of patches was similar to that for the floral traits, using analyses that made use of reference populations. To do so we averaged Q for each patch and then calculated the Spearman's rank correlations between mean Q and the mean CDF score based on floral morphology. To see if the genetic pattern of patches correlated with physical distance, so that differences between the two sites could be explained simply by differences in the distribution of suitable habitat, we also calculated the correlation between mean Q and distance from the end of the zone of sympatry (as shown in Figure 1). Correlations were calculated separately for each zone of sympatry.

As a way of quantifying the genetic structure at each site, we also calculated FST between species using Bayesian maximum likelihood Markov chain Monte Carlo estimation in the program HICKORY (Holsinger et al., 2002). This method estimates population structure using dominant markers, without assuming Hardy–Weinberg's equilibrium to calculate heterozygote frequencies. This analysis required the identification of populations at each site as either I. aggregata or I. tenuituba, and quantified the amount of genetic structure between the species. At GR, all populations were identified as one or the other parent species (I. tenuituba n=7; I. aggregata n=5); at PG, populations A–C were identified as I. tenuituba, and K and L as I. aggregata (Figure 1). Four models were tested: estimation of FST and FIS, estimation of FST assuming FIS=0, estimation of FIS assuming FST=0 and estimation of FST disregarding FIS altogether. A Deviance Information Criterion, which is analogous to Akaike's Information Criterion, was used to choose among models. The hypothesis of greater genetic structure at the bimodal site (GR) was tested by calculating the 95% CI for FST at GR–FST at PG using the results of the Markov chain procedure to generate a null distribution for that parameter.

We compared estimates of genetic structure based on neutral markers (FST from RAPD data) to analogous estimates based on floral morphology. We estimated the structure from floral morphology using QST derived from variance components (McKay and Latta, 2002). A value of QST for a quantitative trait is the proportion of genetic variance that lies at the among-population rather than the within-population level. To obtain estimates of these values, variance components were calculated within each site using the five original morphological traits and employing GLM methods in SPSS, with species as factor. Although our estimates are derived from partitioning phenotypic rather than additive genetic variation, the differences between species in all of the floral traits are retained in common gardens, indicating a strong genetic basis to the differences (Campbell and Waser, 2007). The common garden data for site PG (unavailable for site GR) give rise to similar values of QST to those obtained here for that site. Individuals were assigned to species based on morphology at GR, and assigned to species or ‘hybrid’ based on spatial position at PG.

Results

RAPD markers

Altogether, 242 RAPD primers were screened initially, and 17 used for full screening. Those 17 primers produced 48 polymorphic bands, none of which was species-specific at either site. Most markers were present at similar frequencies at both sites (see Supplementary materials for primer information and marker frequencies). Inheritance was tested for 19 of the 48 markers, none of which deviated significantly from expectations under a model of dominant, Mendelian inheritance (P-values for χ2 tests with d.f.=1 ranged from 0.058 to 1.000). Furthermore, none of the 19 markers showed uniparental inheritance. We could not test inheritance for the other 29 markers because we did not have parent plants that had different genotypes at those loci. We have no reason to assume that any errors associated with RAPD would have been more prevalent at one site than the other, so the comparison of genetic structure among the two sites remains valid even if our genetic data give too conservative an indication of genetic divergence among these species. We note that our RAPD data do generally agree with AFLP data obtained from PG, which we further discuss below.

Morphological structure

Maximum likelihood clustering according to floral morphology produced two clusters at GR and one cluster at PG (Table 1). At GR, the best-fit model was ellipsoidal with two clusters of equal variance, separated clearly according to corolla width (Figure 2a). Those clusters corresponded well to the two species; only four individuals out of 117 total had a probability less than 0.9 of belonging to one or the other cluster. The average uncertainty for classifying individual plants into these two clusters was 0.017, with a median uncertainty of zero. At PG, the best-fit model was ellipsoidal, multivariate normal with one cluster (Figure 2b). All plants were classified in this single cluster with probability one. These clustering results indicated strong bimodality at GR, and unimodality at PG, for floral morphology.

Table 1 Bayesian Information Criterion scores for six multidimensional clustering models tested against empirical data on floral morphology of Ipomopsis aggregata and I. tenuituba from two contact sites in western Colorado, USA
Full size table
Figure 2
figure 2

Representative bivariate plots showing morphometric structures at two Ipomopsis aggregata—I. tenuituba zones of sympatry in Colorado, USA. Morphological data were analyzed with maximum likelihood clustering methods. Plots of normalized corolla length and width show two clusters at Grizzly Ridge, Montrose Co (GR; a) and 1 cluster at Poverty Gulch, Gunnison Co (PG; b); elipses indicate centroids and 95% CI for clusters. Altogether, five morphological traits were used in the analysis; corolla length and width showed the clearest patterns of structure.

Full size image

Genetic structure

The results of clustering test individuals using the program STRUCTURE were quite different for the two contact sites. At GR, K=2 clusters had a posterior probability between 0.9999–1, with all of the I. aggregata patches (based on floral morphology) having an average Q lower than 0.5 and all of the I. tenuituba patches having an average Q higher than 0.5. At PG, K=1 cluster had a posterior probability of 0.9999–1, suggesting an absence of the detectable structure at that site. Thus, this method found two clusters at site GR and one cluster for the hybrid zone itself at site PG.

Incorporating prior information from the reference populations improved the genetic resolution and led to clear differences in the frequency distribution of Q between sites (contingency analysis, G=7.784, d.f.=10, P<0.0001). The resulting frequency distribution of Q was strongly bimodal for site GR (dip=0.0854, n=210, P<0.0001), with few individuals having Q near 0.5 (Figure 3a) and not bimodal for site PG (dip=0.0256, n=215, P>0.4), where many more individuals had intermediate Q values (Figure 3b).

Figure 3
figure 3

Frequency distributions of membership in the cluster defined by the reference population for Ipomopsis tenuituba (Q) at Grizzly Ridge (a) and Poverty Gulch (b). Dark stippled bars: plants in patches classified as I. aggregata based on floral morphology; Light stippled bars: plants in patches classified as hybrid based on floral morphology; Open bars: plants in patches classified as I. tenuituba based on floral morphology.

Full size image

At both sites, the average values for the ancestry value of Q (indicating similarity to I. tenuituba rather than I. aggregata) was strongly correlated with the average value for the CDF based on floral morphology (Spearman's rank r=0.853, N=12 patches, P<0.001 at GR and r=0.812, N=12 patches, P<0.01 at PG). At site PG there was a strong spatial cline in Q, with ancestry in cluster 1 changing gradually from the bottom to the top of the hybrid zone (Figure 4). The correlation of average Q for a patch with distance from the end of the zone of sympatry was nearly perfect (r=−0.926, P<0.0001). The high elevation I. tenuituba patches had high values for Q, whereas the low elevation I. aggregata patches had low Q. Intermediate values for Q were dominated by individuals in the hybrid patches (Figure 3b). In contrast, at site GR there was no detectable correlation of mean Q with distance (r=−0.525, P>0.05), indicating that genetic similarity to I. tenuituba was not closely related to the physical spatial structure at that site.

Figure 4
figure 4

Results of structure using reference populations. Cluster membership Q is given for test individuals at site GR (a) arranged according to classification into the two species based on morphology (I. aggregata vs I. tenuituba) rather than by spatial position. Test individuals at site PG (b) are arranged from the lowest elevation patch to the highest elevation patch.

Full size image

Estimates of genetic structure at the two sites showed weak, but consistent, agreement with the indications of clustering analyses. Bayesian maximum likelihood estimates of FST between species were low at both sites (GR=0.044; PG=0.024), but the 95% CI of FSTGR–FSTPG barely included zero (−0.0022–0.0431), indicating marginal significance for a larger FST at GR than at PG. The full model (estimating FST and FIS) and the FIS=0 model were indistinguishable based on Deviance Information Criterion values, and both were favored over the FST=0 and FIS-free models. The values reported are estimates from the full model. QST for individual floral traits ranged from 0.3 to 0.9, an order of magnitude higher than any of the FST values, and similar to values of QST based on the same traits when the species were raised in common gardens at PG (range, 0.4–1.0; calculated from data in Campbell and Waser, 2007).

Discussion

We found that two zones of sympatry between the same pair of Ipomopsis species differ markedly in frequency of hybridization based on neutral genetic markers. Both RAPD markers and floral morphological traits pointed to contrasting patterns of modality at the two zones of sympatry. Both types of marker indicated a bimodal structure at Grizzly Ridge, where the frequency of apparently hybrid individuals is very low based on corolla color observations and a unimodal distribution at Poverty Gulch where such clearly hybrid individuals are abundant. This suggests that corolla morphology and neutral genetic markers produce similar indications of the rate of hybridization in Ipomopsis aggregata—I. tenuituba zones of sympatry. The contrasting modality at these sites further suggests that hybridization rates vary among zones of sympatry within this species pair. Although our methods did not allow precise identification of the hybrid status (for example, F1 vs F2) of individual plants, they did establish geographic differences in the kind of structure detectable using maximum likelihood clustering methods with no a priori assumptions of structure. These results also agree with previous measures of modality at these sites using hybrid indices calculated from canonical discriminant analysis of floral morphology, which requires assumptions regarding number of clusters and characteristic species morphology (Figure 1—large histograms). Differences in modality between GR and PG were much weaker for neutral genetic markers based on FST than for floral morphology, which is consistent with an interpretation that floral morphology is under selection by pollinators, as has been shown in a previous study (Campbell et al., 1997). Other studies of geographical variation in hybrid frequency in plants have found differences in clinal width (Heywood, 1986) or sharpness of transition between species (Watano et al., 2004), but Buerkle and Rieseberg (2001) found no variation among sites in introgression of particular chromosome blocks in Helianthus.

Introgression of neutral markers in the apparent absence of introgression of morphological traits has been reported elsewhere (Hodges and Arnold, 1994; Vines et al., 2003), and has been interpreted as evidence of diversifying selection for two distinct phenotypes (Barton and Hewitt, 1985). Studies of pollinator preference at PG indicate that at least one of the morphological traits, corolla width, may be under disruptive selection during rare years when not only hummingbirds but also hawkmoths are common pollinators (Campbell et al., 1997). Introgression of neutral traits across a species boundary maintained by selection still requires hybridization and backcrossing, and thus the formation of at least one generation of individuals heterozygous at loci that are under selection. Neutral introgression without morphologically intermediate individuals would be possible if F1s and backcrosses were not necessarily intermediate for selected phenotypic traits, as does happen in many genera (Rieseberg and Ellstrand, 1993). One possible explanation for our results is that the two species have diverged at relatively few loci that control floral morphology and key vegetative and physiological traits (Campbell, 2004), but remain largely homogeneous throughout the rest of the genome.

Variation in modality among zones of sympatry could, in theory, arise in at least three non-mutually exclusive ways. First, it might reflect a simple difference in spatial arrangement of suitable habitat for the two species, such that patches of the two species are closer together at a zone of sympatry that has more hybrids. This explanation can be ruled out, as patches of the two species were actually closer together at Grizzly Ridge, where hybrids were far scarcer. Furthermore, there was no correlation of genotype (assessed by the Q measure of ancestry) with physical distance at that site.

Second, the variation between zones of sympatry could reflect different stages in speciation, due, for example, either to chance differences in appearance of mutations during ongoing primary divergence, or to different timings of secondary contact following allopatric divergence. The existence of zones of sympatry exhibiting a range of reproductive isolation from complete to non-existent has been used as evidence for parapatric speciation (Jiggins and Mallett, 2000), with recent or ongoing in situ evolution of reproductive isolation at each site. This conclusion has been challenged on the grounds that the same pattern could result from secondary contact between many pairs of species that evolved various degrees of reproductive isolation in allopatry (Coyne and Orr, 2004). Variation among zones of sympatry of a single species pair, such as is reported here, controls for among-taxon differences in allopatric divergence unless we allow for different amounts of allopatric divergence among populations of each species, which cannot be ruled out based on our results. It is unclear whether these zones of sympatry are the result of secondary contact or ongoing primary divergence. A broad survey of Ipomopsis zones of sympatry found the presence of site-specific—but absence of species-specific—chloroplast markers, suggesting that these species might be undergoing parapatric speciation. However, those data were difficult to interpret because of a lack of informative nuclear markers (Wolf et al., 1997). A more recent AFLP analysis of PG has produced a single species-specific nuclear marker and several species-specific chloroplast markers for that site (Wu and Campbell, 2005). Should these markers prove to be diagnostic at other Ipomopsis zones as well, the evidence would favor a single speciation event and weigh against the hypothesis that I. aggregata and I. tenuituba have diverged in allopatry or parapatry to different degrees in different places. Wu and Campbell (2005) found a disagreement between distributions of the chloroplast and nuclear markers at PG that suggested I. aggregata nuclear genes might be advancing into populations that have I. tenuituba cytoplasmic genotypes. Our finding of less overall population structure at PG at the nuclear genome level is consistent with that hypothesis.

A third interpretation of our results, and the one with the most evidence at present, is that exogenous reproductive isolation between I. aggregata and I. tenuituba is stronger at GR than at PG because of differences in ecological conditions. This interpretation has considerable support from observations of pollinator behavior. Although hummingbirds are present at both sites, the longer-tongued hawkmoth pollinators were far more common at GR than at PG during the period of study (Aldridge and Campbell, 2007). The reasons for this difference are unknown, and could reflect either consistent historical differences in climatic conditions or a recent anthropogenic change. In any case, the consistent presence of large numbers of both pollinator types at GR could allow for disruptive selection on floral traits as envisioned by Grant and Grant (1965), and thereby explain the lower frequency of hybrids. Indeed, at GR the pollinators exhibit strong floral preferences and rarely fly between plants of the two species even when these are presented intermixed in experimental arrays. Such strong preferences of two pollinators differing in tongue length are predicted by simple optimal foraging models (Rodriguez-Girones and Santamaria, 2006). At PG, in contrast, interspecific flights are more common, perhaps because hawkmoths there forage under higher light intensities when colors of both flower species may be more equally visible (Aldridge and Campbell, 2007). Other factors, such as vegetation structure (Watano et al., 2004) and pollen competition (Williams et al., 2001), have also been suggested to vary among zones. We found no evidence for conspecific pollen advantage at either site, but an asymmetrical seed siring advantage for I. aggregata at GR; however, this advantage alone was insufficient to account for the low frequency of hybrids at GR (Aldridge and Campbell, 2006). Direct comparison of reproductive isolating mechanisms at the two sites thus suggests that modality is strongly influenced by pre-mating reproductive isolation in this species pair, and that local ecological conditions, such as pollinator abundance or behavior or the spatial structure of parent populations (mosaic vs cline), affect the strength of reproductive isolating mechanisms. In such a case, where local conditions affect reproductive isolation, the evolutionary outcome of sympatry between the same two species could be different in different places. Future studies of other systems should investigate potential ecological mechanisms involved in producing such geographical differences in hybridization.