# The Role of Environmental Heterogeneity in Maintaining Reproductive Isolation between HybridizingPasserina (Aves: Cardinalidae) Buntings

Authors: Matthew D. Carling; Henri A. Thomassen
Journal: International Journal of Ecology
 (2012)

---

## Abstract

Hybrid zones are useful systems in which to investigate processes important in creating and maintaining biological diversity. As they are often located in ecotones, patterns of environmental heterogeneity may influence hybridization, and may also influence the maintenance of reproductive isolation between hybridizing species. Focusing on the hybrid zone betweenPasserina amoena (Lazuli Bunting) and Passerina cyanea (Indigo Bunting), located in the eastern Rocky Mountain/western Great Plains ecotone, we examined the relationship between population-pairwise differences in the proportion of hybrids and environmental variation. Models including environmental variables explained more of the variation in hybridization rates across the ecotone than did models that only included the geographic distance between sampling localities as predictor variables (63.9% and 58.9% versus 38.8% and 39.9%, depending on how hybridization was quantified). In the models including environmental variables, the amount of rainfall during the warmest quarter had the greatest explanatory power, consistent with a hypothesis that P. cyanea is better adapted to the mesic environments of eastern North America and P. amoena is better adapted to the xeric habitats of western North America. These results suggest that continued reproductive isolation between these species is mediated, at least partially, by differential adaptations to local environmental conditions.


---

## Body

## 1. Introduction

Hybrid zones formed through the interbreeding of divergent lineages as a result of secondary contact are powerful models for investigating the forces important in generating and maintaining biodiversity [1–5]. These “natural laboratories” [1] provide an opportunity to explore patterns of gene flow and introgression between closely related taxa, which can aid in the identification of genes contributing to reproductive isolation [6, 7]. In addition, naturally occurring hybrid zones allow the investigation of how and to what extent ecological differences between the hybridizing lineages contribute to speciation [8, 9].Divergent natural selection between different habitats can be an important force in promoting speciation through its potential to cause habitat isolation between closely related taxa [10–12]. When traits subject to divergent natural selection are associated or correlated with traits that may influence reproduction, ecological speciation may result [12]. Naturally occurring hybrid zones are often located in ecotones produced by sharp transitions between different environments in which the hybridizing species were subjected to divergent natural selection [13, 14]. Such ecotones are often characterized by areas of marginal habitat that may not be ideal for the parental species found on either side of the environmental gradient and therefore provide an opportunity to examine the impact of habitat on hybridization and the evolution of reproductive isolation [8, 11]. In this manner, hybrid zones offer a chance to explore the environmental conditions that were potentially important to the initial divergence of the hybridizing species, providing insight into the importance of ecological speciation, mediated through divergent natural selection.Spatial environmental heterogeneity can change over time, with the potential to greatly alter the structure of hybrid zones, and influence speciation and biodiversity dynamics [8]. If hybridizing species are differentially adapted to the environments they primarily inhabit, changes in those environments and the intervening ecotone may result in one of two general outcomes: (1) reduction in species diversity—one species may swamp the other, causing extinction of the overwhelmed species, or the species may collapse into a hybrid swarm wherein both parental species in essence become extinct, but a new “hybrid” species is formed—or (2) reduction or cessation of hybridization—if hybridization only occurs in the “hybrid” habitat found in an ecotone and such habitats are altered or eliminated. Alternatively, environmental stability may promote the stability of a hybrid zone, regardless of whether the zone is maintained through endogenous selection against hybrids balanced by continued migration into the zone (i.e., a tension zone) or exogenous selection creating fitness differences along the environmental gradient found in the ecotone.In this paper, we quantify the relationship between patterns of environmental heterogeneity and patterns of hybridization across an ecotone, using the hybrid zone formed between Lazuli (Passerina amoena) and Indigo (Passerina cyanea) buntings as a model. We explore whether differences between populations in the proportion of hybrid individuals are better predicted by the variation found along the ecological gradient experienced across the hybrid zone than by simple geographic distance between populations. In doing so, we investigate the potential role of ecological differences between P. amoena and P. cyanea in maintaining reproductive isolation between these species.Passerina amoena and P. cyanea are closely related passerines [15] that form a hybrid zone in the ecotone located between the eastern Rocky Mountains and the western Great Plains of North America (Figure 1; [16–18]). Behavioral studies, carried out in both the lab and the field, demonstrate strong patterns of assortative mating and show that, when choosing mates, females of both species place greater emphasis on male plumage patterns than on song characteristics [19–21]. Field-based research also suggests that hybrid individuals suffer fitness consequences relative to pure individuals. Mated pairs that involve at least one hybrid individual produce fewer nestling and fledglings than do pairs involving only nonhybrid individuals [21].(a) Primary breeding distributions and allopatric sampling localities ofPasserina amoena and Passerina cyanea. Area outlined by the black box is enlarged in (b). (b) Sympatric and parapatric sampling localities analyzed in this study. Background shading represents elevation and white box indicates approximate location of the sympatric contact zone. See Table 1 for more detailed information on sampling localities. Digital maps for each species [31] were downloaded from NatureServe [32] and modified.
(a)(b)More recently, coalescent analyses of the evolutionary history ofP. amoena and P. cyanea suggest that speciation occurred through either parapatric divergence or repeated cycles of allopatric divergence with gene flow during periods of secondary contact [22]. Regardless of the spatial context of divergence, the current hybrid zone originated through secondary contact and is probably no older than ~6500 years [7, 22, 23]. Patterns of genetic introgression across the hybrid zone, based on allelic data from two mitochondrial genes, four nuclear autosomal loci, and two sex-linked loci, indicate reduced introgression of mitochondrial DNA and sex-linked loci, relative to the pattern seen for nuclear autosomal loci [23], a finding consistent with Haldane’s rule [24].Additional work has shown that introgression patterns of sex-linked loci can vary substantially [7]; maximum likelihood estimates of 10 locus-specific cline widths differed by ~200-fold, identifying a candidate genomic region for reproductive isolation. Suggestively, the sex-linked locus with the narrowest cline width is an intron in the very-low-density lipoprotein receptor gene (VLDLR); mutations in this gene are known to reduce egg laying capabilities in chickens [25] and it is possible that divergent P. cyanea and P. amoena VLDLR alleles do not function well in a heterospecific genomic background, which could cause hybrid females to lay fewer eggs. The hypothesis that a Dobzhansky-Muller incompatibility [26–28] causes a reduction in egg laying in hybrid females is consistent with the field research findings of Baker and Boylan [21] described above.Using a species distribution modeling approach, Swenson [29] suggested that the ecotonal environment between the eastern Rocky Mountains and western Great Plains is important in determining the geographic location of the Passerina bunting hybrid zone, with temperature differences having the biggest impact. This may partially explain the observed shift in location and width of the hybrid zone that has occurred over the past 40–45 years [30]. However, Swenson’s [29] analysis utilized only the geographic location of hybrid individuals; here we quantify the proportion of hybrid individuals within a population to elucidate the influence of environmental heterogeneity on hybridization dynamics.


## 2. Materials and Methods

### 2.1. Genetic Data and Analyses

Population samples ofP. amoena, P. cyanea, and P. amoena × P. cyanea hybrids were collected from 21 sympatric and parapatric localities across the contact zone during May, June, and July 2004–2007 (Figure 1, Table 1). All collected individuals were prepared as voucher specimens and deposited in the Louisiana State University Museum of Vertebrates. Additional samples from allopatric populations (P. amoena: WA, ID; P. cyanea: MN, IL, MI) were also obtained (Figure 1, Table 1).Table 1
Sampling localities, sample sizes, and the proportion of hybrid individuals for populations analyzed in this study.
CodeLocalityLatitudeLongitudeNProportion of hybrids (80–20)Proportion of hybrids (90–10)1WA47.8−118.929——2ID44.1−116.23——3MT: Custer National Forest no. 145.1−108.540.000.004WY: Bighorn National Forest no. 144.8−107.380.000.005WY: Bighorn National Forest no. 244.6−107.720.000.006CO: Roosevelt National Forest40.7−105.270.140.147MT: Custer National Forest no. 245.7−10650.000.008WY: Medicine Bow National Forest42.4−105.330.000.009ND: Little Missouri National Grassland46.8−103.5130.000.0010SD: Custer National Forest45.3−103.2100.000.0011WY: Sand Creek*44.5−104.1140.360.5012NE: White River*42.6−103.580.380.5013SD: Black Hills National Forest43.7−103.850.200.2014NE: Ponderosa State Wildlife Management Area42.6−103.370.571.0015SD: The Nature Conservancy Whitney Preserve43.3−103.630.000.0016NE: Nebraska National Forest42.8−102.980.251.0017SD: Ft. Meade National Recreation Area44.4−103.540.500.5018NE: Nenzel*42.8−101.140.250.2519NE: The Nature Conservancy42.8−100180.000.00Niobrara Valley Preserve20SD: Carpenter Game Production Area43.7−99.5130.080.0821ND: The Nature Conservancy46.3−97.340.000.00Pigeon Point Preserve22NE: Wiseman State Wildlife Area42.8−97.1100.000.0023SD: Newton Hills Game Production Area43.2−96.6130.000.0824MN44.9−93.710——25IL40.5−88.910——26MI42.3−83.77——*Access to these localities generously provided by private landowners.From each individual, genomic DNA was extracted from pectoral muscle and used to generate sequence data from a suite of 14 loci (see Supplementary Table in Supplementary Material available online at doi: 10.1155/2011/295463 Table  S1) as previously described [7, 22]. We resolved individual haplotypes probabilistically using PHASE [33, 34] and then identified the largest independently segregating block of sequence for all loci as in Carling and Brumfield [7, 22]. Individual haplotypes with frequencies ≥0.80 in the “allopatric” P. amoena population (WA) were classified as P. amoena haplotypes, and haplotypes with frequencies ≥0.80 in the “allopatric” P. cyanea populations (MN, IL, MI) were designated as P. cyaneahaplotypes [7, 22]. We then generated multilocus genotypes for each individual by coding each allele at each locus as either P. amoenaor P. cyanea.Using these multilocus genotypes we quantified the genetic ancestry of each individual using STRUCTURE [35–38]. In the analyses, the allopatric populations (P. amoena: WA, ID; P. cyanea: MN, IL, MI) were coded as “learning populations” (USEPOPINFO model). For all other individuals, we estimated the proportion of their ancestry attributed to the two clusters defined by the “learning populations.” We used the default allele frequency model and ran the analysis three times, each time with a burnin length of 100,000 steps followed by a postburnin run of 1,000,000 steps. We then calculated the mean ancestry values across the three runs for each individual from the sympatric and parapatric populations. Individuals were classified as a “hybrid” according to two different criteria: (1) “80–20” wherein an individual was classified as a hybrid if the proportion of its genetic ancestry assigned to cluster 1 (P. amoena) was between 0.80 and 0.20 and (2) “90–10” wherein an individual was classified as a hybrid if the proportion of its genetic ancestry assigned to cluster 1 was between 0.90 and 0.10. In both scenarios, individuals whose proportion of genetic ancestry assigned to cluster 1 was greater than 0.80 (criterion 1) or 0.90 (criterion 2) were classified as pure P. amoena individuals and all other individuals were classified as pure P. cyaneaindividuals.Within each sympatric and parapatric population, we calculated the proportion of hybrid individuals, which were then used to generate a distance matrix of the relative similarity in hybridization among populations. The “hybridization distance” between two populations was calculated as(1)|(PL1+0.5PH1)-(PL2+0.5PH2)|,
where PL1 is the proportion of P. amoena individuals in population 1, PH1 is the proportion of hybrids in population 1, PL2 is the proportion of P. amoena individuals in population 2, and PH2 is the proportion of hybrids in population 2. This hybridization distance was calculated for every pair of populations, except the allopatric ones (WA, ID, MN, IL, MI) and the resultant matrix (Table 2) formed the basis of our generalized dissimilarity modeling (see below). As an alternative measure of “hybridization distance” between two populations, we also calculated the absolute difference between mean population Q-scores, as estimated by STRUCTURE [35–38].Table 2
Matrix of hybridization distances(|(PL1+0.5PH1)-(PL2+0.5PH2)|; see text for details) between all pairs of sympatric and allopatric populations. Above the diagonal, hybrids were identified using the “80–20” criterion; below the diagonal, hybrids were identified using the “90–10” criterion.
MT: Custer NF no.1WY: Bighorn NF no.1WY: Bighorn NF no.2CO: Roosevelt NFMT: Custer NF no.2WY: Medicine Bow NFND: Little Missouri NGSD: Custer NFWY: Sand CreekNE: White RiverSD: Black Hills NFNE: Ponderosa SWMASD: TNC Whitney PreserveNE: Nebraska NFSD: Ft. Meade NRANE: NenzelNE: TNC Niobrara Valley PreserveSD: Carpenter GPAND: TNC Pigeon Point PreserveNE: Wiseman SWASD: Newton Hills GPAMT: Custer NF no. 1—0.000.000.070.000.000.000.000.320.190.100.290.000.130.500.881.000.961.001.001.00WY: Bighorn NF no. 10.00—0.000.070.000.000.000.000.320.190.100.290.000.130.500.881.000.961.001.001.00WY: Bighorn NF no. 20.000.00—0.000.000.000.000.000.320.190.100.290.000.130.500.881.000.961.001.001.00CO: Roosevelt NF0.070.070.07—0.070.070.070.070.250.120.030.210.070.050.430.800.930.890.930.930.93MT: Custer NF no. 20.000.000.000.07—0.000.000.000.320.190.100.290.000.130.500.881.000.961.001.001.00WY: Medicine Bow NF0.000.000.000.070.00—0.000.000.320.190.100.290.000.130.500.881.000.961.001.001.00ND: Little Missouri NG0.000.000.000.070.000.00—0.000.320.190.100.290.000.130.500.881.000.961.001.001.00SD: Custer NF0.000.000.000.070.000.000.00—0.320.190.100.290.000.130.500.881.000.961.001.001.00WY: Sand Creek0.320.320.320.250.320.320.320.32—0.130.220.040.320.200.180.550.680.640.680.680.68NE: White River0.250.250.250.180.250.250.250.250.07—0.090.100.190.060.310.690.810.770.810.810.81SD: Black Hills NF0.100.100.100.030.100.100.100.100.220.15—0.190.100.030.400.780.900.860.900.900.90NE: Ponderosa SWMA0.500.500.500.430.500.500.500.500.180.250.40—0.290.160.210.590.710.680.710.710.71SD: TNC Whitney Preserve0.000.000.000.070.000.000.000.000.320.250.100.50—0.130.500.881.000.961.001.001.00NE: Nebraska NF0.500.500.500.430.500.500.500.500.180.250.400.000.50—0.380.750.880.840.880.880.88SD: Ft. Meade NRA0.500.500.500.430.500.500.500.500.180.250.400.000.500.00—0.380.500.460.500.500.50NE: Nenzel0.880.880.880.800.880.880.880.880.550.630.780.380.880.380.38—0.130.090.130.130.13NE: TNC Niobrara Valley Preserve1.001.001.000.931.001.001.001.000.680.750.900.501.000.500.500.13—0.040.000.000.00SD: Carpenter GPA0.960.960.960.890.960.960.960.960.640.710.860.460.960.460.460.090.04—0.040.040.04ND: TNC Pigeon Point            Preserve1.001.001.000.931.001.001.001.000.680.750.900.501.000.500.500.130.000.04—0.000.00NE: Wiseman SWA1.001.001.000.931.001.001.001.000.680.750.900.501.000.500.500.130.000.040.00—0.00SD: Newton Hills GPA0.960.960.960.890.960.960.960.960.640.710.860.460.960.460.460.090.040.000.040.04—


### 2.2. Spatial Data and Analyses

#### 2.2.1. Environmental Data

We used a set of moderately high-resolution climate and satellite remote sensing variables to characterize the habitat within and on both sides of the hybrid zone (Table3). Climate data were obtained from the WorldClim database [39], which are spatially explicit estimates of annual means, seasonal extremes, and degrees of seasonality in temperature and precipitation based on a 50-year climatology (1950–2000). The climate variables that were included in our analyses are annual mean temperature (Bio1); mean diurnal temperature range (Bio2), which is the difference in daily maximum and minimum temperatures, averaged across the year; mean temperature of the warmest quarter (Bio10);mean temperature of the coldest quarter (Bio11);total annual precipitation (Bio12); precipitation seasonality (Bio15), which is the coefficient of variation in monthly rainfall across the year;total precipitation of the wettest quarter (Bio17); total precipitation of the warmest quarter (Bio18).Table 3
List of environmental variables used in our analyses.
Data RecordInstrumentEcological attributesVariables derivedNormalized Difference Vegetation Index (NDVI)Satellite-MODISVegetation densityNDVI meanVegetation Continuous Field (VCF)*§Satellite-MODISForest cover + heterogeneityTree coverScatterometer-Backscatter*†Satellite-QSCATSurface moisture + roughness (forest structure), seasonalityQSCAT  augDEM*SpaceShuttle SRTMTopography + ruggednessSRTMWorldClim¶Station-NetworkBioclimatic variablesBio1, Bio2, Bio10, Bio11, Bio12, Bio15, Bio17, Bio18*Data at native resolutions smaller or larger than 1 km have been aggregated to 1 km.†QSCAT is based on data from 2001.§Based on MODIS data from 2001 [41].¶WorldClim data are based on monthly climatologies from 1950–2000 [39]. The bioclimatic variables are annual mean temperature (Bio1), mean diurnal temperature range (Bio2), mean temperature of the warmest quarter (Bio10), mean temperature of the coldest quarter (Bio11), annual precipitation (Bio12), precipitation seasonality (Bio15), precipitation of the wettest quarter (Bio17), and precipitation of the warmest quarter (Bio18).In addition to these ground-based measurements of climate, we used satellite remote sensing data from both passive optical sensors (MODIS;https://lpdaac.usgs.gov/lpdaac/products/modis_overview) and active radar scatterometers (QuickScat; http://www.scp.byu.edu/data/Quikscat/SIRv2/qush/World_regions.htm) to infer a variety of ecological characteristics of the land surface. From the MODIS archive, we used the average monthly Normalized Difference Vegetation Index (NDVI) to infer vegetation density [40]. NDVI is sometimes also referred to as a measure of “greenness” and is high over forested areas and low over areas with sparse vegetation. In addition, we used the vegetation continuous field [41] product as a measure of the percentage of tree cover. In contrast to NDVI, percent tree cover specifically measures the land surface area covered by trees, disregarding other types of vegetation. From QuickScat (QSCAT), we obtained raw backscatter measurements for the month of August that capture attributes related to surface moisture [40], and, from the Shuttle Radar Topography Mission (SRTM), we acquired elevation data.Time series of the remote sensing data sources were acquired from 2001 for QSCAT and tree cover and means over the period 2000–2004 for NDVI. Variables with native resolutions higher (e.g., SRTM: 30 m) or lower (e.g., QSCAT: 2.25 km) than 1 km were reaggregated to a 1 km grid cell resolution. This resolution is often used in ecological niche modeling at the regional or subcontinental scales and balances resolutions from climate data, which often have coarser native resolutions, and remote sensing data, with higher native resolutions. To improve interpretation, we checked for covariance among variables and only included those with substantial unique variance. Various criteria were used to decide which layers of correlated pairs (with Pearson's correlations on the order of 0.9 or larger) were retained for further analysis. These included keeping layers that are more commonly used in distribution modeling (WorldClim) or that exhibit larger contrast/variance over the study area (QSCAT) as well as having best data quality (NDVI). Pearson’s cross correlations of the used environmental variables are shown in supplementary Table S2.


#### 2.2.2. Generalized Dissimilarity Modeling

To analyze the relative contribution of geographic distance and environmental variables to explaining the similarity among sampling locations in the percentage of hybrids, we used generalized dissimilarity modeling (GDM, [42]). GDM is a matrix regression technique that predicts biotic dissimilarity (turnover) between sites based upon environmental heterogeneity and geographic distance. One advantage of GDM over other modeling methodologies is that it can fit nonlinear relationships between predictor and response variables through the use of I-spline basis functions [42]. Another advantage of GDM is that, because it uses pairwise comparisons between sites, it can explicitly take into account the influence of geographic distance on explaining biological variation—a particularly important feature in our study. The relative importance of predictor variables in a GDM can be assessed by means of response curves. To further evaluate the extent to which geographic distance is potentially correlated with environmental differences, we ran independent tests with the following sets of predictor variables: (1) environmental variables and geographic distance; (2) only geographic distance; (3) only environmental variables.GDM is a two-step method: first, dissimilarities of a set of predictor variables are fitted to the genetic dissimilarities (the response variables). The contributions of predictor variables to explaining the observed response variation are tested by permutations, and only those variables that are significant are retained in the final model. These procedures result in a function that describes the relationship between environmental and response variables. Second, using the function resulting from the first step, a spatial prediction is made of the response variable patterns. For visualization purposes, classes of similar response are color coded, where larger color differences between two localities represent larger differences in the proportion of hybrid individuals. In order to assess the significance of the level of variation that was explained by our models, we ran additional models in which the environmental layers were substituted by layers with random values for each grid cell. The resulting percentage of variation explained was compared to that of the full model. We considered the performance of the full model not significant if it explained an equal amount or less of the total variation than a model with random environmental variables.


## 2.1. Genetic Data and Analyses

Population samples ofP. amoena, P. cyanea, and P. amoena × P. cyanea hybrids were collected from 21 sympatric and parapatric localities across the contact zone during May, June, and July 2004–2007 (Figure 1, Table 1). All collected individuals were prepared as voucher specimens and deposited in the Louisiana State University Museum of Vertebrates. Additional samples from allopatric populations (P. amoena: WA, ID; P. cyanea: MN, IL, MI) were also obtained (Figure 1, Table 1).Table 1
Sampling localities, sample sizes, and the proportion of hybrid individuals for populations analyzed in this study.
CodeLocalityLatitudeLongitudeNProportion of hybrids (80–20)Proportion of hybrids (90–10)1WA47.8−118.929——2ID44.1−116.23——3MT: Custer National Forest no. 145.1−108.540.000.004WY: Bighorn National Forest no. 144.8−107.380.000.005WY: Bighorn National Forest no. 244.6−107.720.000.006CO: Roosevelt National Forest40.7−105.270.140.147MT: Custer National Forest no. 245.7−10650.000.008WY: Medicine Bow National Forest42.4−105.330.000.009ND: Little Missouri National Grassland46.8−103.5130.000.0010SD: Custer National Forest45.3−103.2100.000.0011WY: Sand Creek*44.5−104.1140.360.5012NE: White River*42.6−103.580.380.5013SD: Black Hills National Forest43.7−103.850.200.2014NE: Ponderosa State Wildlife Management Area42.6−103.370.571.0015SD: The Nature Conservancy Whitney Preserve43.3−103.630.000.0016NE: Nebraska National Forest42.8−102.980.251.0017SD: Ft. Meade National Recreation Area44.4−103.540.500.5018NE: Nenzel*42.8−101.140.250.2519NE: The Nature Conservancy42.8−100180.000.00Niobrara Valley Preserve20SD: Carpenter Game Production Area43.7−99.5130.080.0821ND: The Nature Conservancy46.3−97.340.000.00Pigeon Point Preserve22NE: Wiseman State Wildlife Area42.8−97.1100.000.0023SD: Newton Hills Game Production Area43.2−96.6130.000.0824MN44.9−93.710——25IL40.5−88.910——26MI42.3−83.77——*Access to these localities generously provided by private landowners.From each individual, genomic DNA was extracted from pectoral muscle and used to generate sequence data from a suite of 14 loci (see Supplementary Table in Supplementary Material available online at doi: 10.1155/2011/295463 Table  S1) as previously described [7, 22]. We resolved individual haplotypes probabilistically using PHASE [33, 34] and then identified the largest independently segregating block of sequence for all loci as in Carling and Brumfield [7, 22]. Individual haplotypes with frequencies ≥0.80 in the “allopatric” P. amoena population (WA) were classified as P. amoena haplotypes, and haplotypes with frequencies ≥0.80 in the “allopatric” P. cyanea populations (MN, IL, MI) were designated as P. cyaneahaplotypes [7, 22]. We then generated multilocus genotypes for each individual by coding each allele at each locus as either P. amoenaor P. cyanea.Using these multilocus genotypes we quantified the genetic ancestry of each individual using STRUCTURE [35–38]. In the analyses, the allopatric populations (P. amoena: WA, ID; P. cyanea: MN, IL, MI) were coded as “learning populations” (USEPOPINFO model). For all other individuals, we estimated the proportion of their ancestry attributed to the two clusters defined by the “learning populations.” We used the default allele frequency model and ran the analysis three times, each time with a burnin length of 100,000 steps followed by a postburnin run of 1,000,000 steps. We then calculated the mean ancestry values across the three runs for each individual from the sympatric and parapatric populations. Individuals were classified as a “hybrid” according to two different criteria: (1) “80–20” wherein an individual was classified as a hybrid if the proportion of its genetic ancestry assigned to cluster 1 (P. amoena) was between 0.80 and 0.20 and (2) “90–10” wherein an individual was classified as a hybrid if the proportion of its genetic ancestry assigned to cluster 1 was between 0.90 and 0.10. In both scenarios, individuals whose proportion of genetic ancestry assigned to cluster 1 was greater than 0.80 (criterion 1) or 0.90 (criterion 2) were classified as pure P. amoena individuals and all other individuals were classified as pure P. cyaneaindividuals.Within each sympatric and parapatric population, we calculated the proportion of hybrid individuals, which were then used to generate a distance matrix of the relative similarity in hybridization among populations. The “hybridization distance” between two populations was calculated as(1)|(PL1+0.5PH1)-(PL2+0.5PH2)|,
where PL1 is the proportion of P. amoena individuals in population 1, PH1 is the proportion of hybrids in population 1, PL2 is the proportion of P. amoena individuals in population 2, and PH2 is the proportion of hybrids in population 2. This hybridization distance was calculated for every pair of populations, except the allopatric ones (WA, ID, MN, IL, MI) and the resultant matrix (Table 2) formed the basis of our generalized dissimilarity modeling (see below). As an alternative measure of “hybridization distance” between two populations, we also calculated the absolute difference between mean population Q-scores, as estimated by STRUCTURE [35–38].Table 2
Matrix of hybridization distances(|(PL1+0.5PH1)-(PL2+0.5PH2)|; see text for details) between all pairs of sympatric and allopatric populations. Above the diagonal, hybrids were identified using the “80–20” criterion; below the diagonal, hybrids were identified using the “90–10” criterion.
MT: Custer NF no.1WY: Bighorn NF no.1WY: Bighorn NF no.2CO: Roosevelt NFMT: Custer NF no.2WY: Medicine Bow NFND: Little Missouri NGSD: Custer NFWY: Sand CreekNE: White RiverSD: Black Hills NFNE: Ponderosa SWMASD: TNC Whitney PreserveNE: Nebraska NFSD: Ft. Meade NRANE: NenzelNE: TNC Niobrara Valley PreserveSD: Carpenter GPAND: TNC Pigeon Point PreserveNE: Wiseman SWASD: Newton Hills GPAMT: Custer NF no. 1—0.000.000.070.000.000.000.000.320.190.100.290.000.130.500.881.000.961.001.001.00WY: Bighorn NF no. 10.00—0.000.070.000.000.000.000.320.190.100.290.000.130.500.881.000.961.001.001.00WY: Bighorn NF no. 20.000.00—0.000.000.000.000.000.320.190.100.290.000.130.500.881.000.961.001.001.00CO: Roosevelt NF0.070.070.07—0.070.070.070.070.250.120.030.210.070.050.430.800.930.890.930.930.93MT: Custer NF no. 20.000.000.000.07—0.000.000.000.320.190.100.290.000.130.500.881.000.961.001.001.00WY: Medicine Bow NF0.000.000.000.070.00—0.000.000.320.190.100.290.000.130.500.881.000.961.001.001.00ND: Little Missouri NG0.000.000.000.070.000.00—0.000.320.190.100.290.000.130.500.881.000.961.001.001.00SD: Custer NF0.000.000.000.070.000.000.00—0.320.190.100.290.000.130.500.881.000.961.001.001.00WY: Sand Creek0.320.320.320.250.320.320.320.32—0.130.220.040.320.200.180.550.680.640.680.680.68NE: White River0.250.250.250.180.250.250.250.250.07—0.090.100.190.060.310.690.810.770.810.810.81SD: Black Hills NF0.100.100.100.030.100.100.100.100.220.15—0.190.100.030.400.780.900.860.900.900.90NE: Ponderosa SWMA0.500.500.500.430.500.500.500.500.180.250.40—0.290.160.210.590.710.680.710.710.71SD: TNC Whitney Preserve0.000.000.000.070.000.000.000.000.320.250.100.50—0.130.500.881.000.961.001.001.00NE: Nebraska NF0.500.500.500.430.500.500.500.500.180.250.400.000.50—0.380.750.880.840.880.880.88SD: Ft. Meade NRA0.500.500.500.430.500.500.500.500.180.250.400.000.500.00—0.380.500.460.500.500.50NE: Nenzel0.880.880.880.800.880.880.880.880.550.630.780.380.880.380.38—0.130.090.130.130.13NE: TNC Niobrara Valley Preserve1.001.001.000.931.001.001.001.000.680.750.900.501.000.500.500.13—0.040.000.000.00SD: Carpenter GPA0.960.960.960.890.960.960.960.960.640.710.860.460.960.460.460.090.04—0.040.040.04ND: TNC Pigeon Point            Preserve1.001.001.000.931.001.001.001.000.680.750.900.501.000.500.500.130.000.04—0.000.00NE: Wiseman SWA1.001.001.000.931.001.001.001.000.680.750.900.501.000.500.500.130.000.040.00—0.00SD: Newton Hills GPA0.960.960.960.890.960.960.960.960.640.710.860.460.960.460.460.090.040.000.040.04—


## 2.2. Spatial Data and Analyses

### 2.2.1. Environmental Data

We used a set of moderately high-resolution climate and satellite remote sensing variables to characterize the habitat within and on both sides of the hybrid zone (Table3). Climate data were obtained from the WorldClim database [39], which are spatially explicit estimates of annual means, seasonal extremes, and degrees of seasonality in temperature and precipitation based on a 50-year climatology (1950–2000). The climate variables that were included in our analyses are annual mean temperature (Bio1); mean diurnal temperature range (Bio2), which is the difference in daily maximum and minimum temperatures, averaged across the year; mean temperature of the warmest quarter (Bio10);mean temperature of the coldest quarter (Bio11);total annual precipitation (Bio12); precipitation seasonality (Bio15), which is the coefficient of variation in monthly rainfall across the year;total precipitation of the wettest quarter (Bio17); total precipitation of the warmest quarter (Bio18).Table 3
List of environmental variables used in our analyses.
Data RecordInstrumentEcological attributesVariables derivedNormalized Difference Vegetation Index (NDVI)Satellite-MODISVegetation densityNDVI meanVegetation Continuous Field (VCF)*§Satellite-MODISForest cover + heterogeneityTree coverScatterometer-Backscatter*†Satellite-QSCATSurface moisture + roughness (forest structure), seasonalityQSCAT  augDEM*SpaceShuttle SRTMTopography + ruggednessSRTMWorldClim¶Station-NetworkBioclimatic variablesBio1, Bio2, Bio10, Bio11, Bio12, Bio15, Bio17, Bio18*Data at native resolutions smaller or larger than 1 km have been aggregated to 1 km.†QSCAT is based on data from 2001.§Based on MODIS data from 2001 [41].¶WorldClim data are based on monthly climatologies from 1950–2000 [39]. The bioclimatic variables are annual mean temperature (Bio1), mean diurnal temperature range (Bio2), mean temperature of the warmest quarter (Bio10), mean temperature of the coldest quarter (Bio11), annual precipitation (Bio12), precipitation seasonality (Bio15), precipitation of the wettest quarter (Bio17), and precipitation of the warmest quarter (Bio18).In addition to these ground-based measurements of climate, we used satellite remote sensing data from both passive optical sensors (MODIS;https://lpdaac.usgs.gov/lpdaac/products/modis_overview) and active radar scatterometers (QuickScat; http://www.scp.byu.edu/data/Quikscat/SIRv2/qush/World_regions.htm) to infer a variety of ecological characteristics of the land surface. From the MODIS archive, we used the average monthly Normalized Difference Vegetation Index (NDVI) to infer vegetation density [40]. NDVI is sometimes also referred to as a measure of “greenness” and is high over forested areas and low over areas with sparse vegetation. In addition, we used the vegetation continuous field [41] product as a measure of the percentage of tree cover. In contrast to NDVI, percent tree cover specifically measures the land surface area covered by trees, disregarding other types of vegetation. From QuickScat (QSCAT), we obtained raw backscatter measurements for the month of August that capture attributes related to surface moisture [40], and, from the Shuttle Radar Topography Mission (SRTM), we acquired elevation data.Time series of the remote sensing data sources were acquired from 2001 for QSCAT and tree cover and means over the period 2000–2004 for NDVI. Variables with native resolutions higher (e.g., SRTM: 30 m) or lower (e.g., QSCAT: 2.25 km) than 1 km were reaggregated to a 1 km grid cell resolution. This resolution is often used in ecological niche modeling at the regional or subcontinental scales and balances resolutions from climate data, which often have coarser native resolutions, and remote sensing data, with higher native resolutions. To improve interpretation, we checked for covariance among variables and only included those with substantial unique variance. Various criteria were used to decide which layers of correlated pairs (with Pearson's correlations on the order of 0.9 or larger) were retained for further analysis. These included keeping layers that are more commonly used in distribution modeling (WorldClim) or that exhibit larger contrast/variance over the study area (QSCAT) as well as having best data quality (NDVI). Pearson’s cross correlations of the used environmental variables are shown in supplementary Table S2.


### 2.2.2. Generalized Dissimilarity Modeling

To analyze the relative contribution of geographic distance and environmental variables to explaining the similarity among sampling locations in the percentage of hybrids, we used generalized dissimilarity modeling (GDM, [42]). GDM is a matrix regression technique that predicts biotic dissimilarity (turnover) between sites based upon environmental heterogeneity and geographic distance. One advantage of GDM over other modeling methodologies is that it can fit nonlinear relationships between predictor and response variables through the use of I-spline basis functions [42]. Another advantage of GDM is that, because it uses pairwise comparisons between sites, it can explicitly take into account the influence of geographic distance on explaining biological variation—a particularly important feature in our study. The relative importance of predictor variables in a GDM can be assessed by means of response curves. To further evaluate the extent to which geographic distance is potentially correlated with environmental differences, we ran independent tests with the following sets of predictor variables: (1) environmental variables and geographic distance; (2) only geographic distance; (3) only environmental variables.GDM is a two-step method: first, dissimilarities of a set of predictor variables are fitted to the genetic dissimilarities (the response variables). The contributions of predictor variables to explaining the observed response variation are tested by permutations, and only those variables that are significant are retained in the final model. These procedures result in a function that describes the relationship between environmental and response variables. Second, using the function resulting from the first step, a spatial prediction is made of the response variable patterns. For visualization purposes, classes of similar response are color coded, where larger color differences between two localities represent larger differences in the proportion of hybrid individuals. In order to assess the significance of the level of variation that was explained by our models, we ran additional models in which the environmental layers were substituted by layers with random values for each grid cell. The resulting percentage of variation explained was compared to that of the full model. We considered the performance of the full model not significant if it explained an equal amount or less of the total variation than a model with random environmental variables.


## 2.2.1. Environmental Data

We used a set of moderately high-resolution climate and satellite remote sensing variables to characterize the habitat within and on both sides of the hybrid zone (Table3). Climate data were obtained from the WorldClim database [39], which are spatially explicit estimates of annual means, seasonal extremes, and degrees of seasonality in temperature and precipitation based on a 50-year climatology (1950–2000). The climate variables that were included in our analyses are annual mean temperature (Bio1); mean diurnal temperature range (Bio2), which is the difference in daily maximum and minimum temperatures, averaged across the year; mean temperature of the warmest quarter (Bio10);mean temperature of the coldest quarter (Bio11);total annual precipitation (Bio12); precipitation seasonality (Bio15), which is the coefficient of variation in monthly rainfall across the year;total precipitation of the wettest quarter (Bio17); total precipitation of the warmest quarter (Bio18).Table 3
List of environmental variables used in our analyses.
Data RecordInstrumentEcological attributesVariables derivedNormalized Difference Vegetation Index (NDVI)Satellite-MODISVegetation densityNDVI meanVegetation Continuous Field (VCF)*§Satellite-MODISForest cover + heterogeneityTree coverScatterometer-Backscatter*†Satellite-QSCATSurface moisture + roughness (forest structure), seasonalityQSCAT  augDEM*SpaceShuttle SRTMTopography + ruggednessSRTMWorldClim¶Station-NetworkBioclimatic variablesBio1, Bio2, Bio10, Bio11, Bio12, Bio15, Bio17, Bio18*Data at native resolutions smaller or larger than 1 km have been aggregated to 1 km.†QSCAT is based on data from 2001.§Based on MODIS data from 2001 [41].¶WorldClim data are based on monthly climatologies from 1950–2000 [39]. The bioclimatic variables are annual mean temperature (Bio1), mean diurnal temperature range (Bio2), mean temperature of the warmest quarter (Bio10), mean temperature of the coldest quarter (Bio11), annual precipitation (Bio12), precipitation seasonality (Bio15), precipitation of the wettest quarter (Bio17), and precipitation of the warmest quarter (Bio18).In addition to these ground-based measurements of climate, we used satellite remote sensing data from both passive optical sensors (MODIS;https://lpdaac.usgs.gov/lpdaac/products/modis_overview) and active radar scatterometers (QuickScat; http://www.scp.byu.edu/data/Quikscat/SIRv2/qush/World_regions.htm) to infer a variety of ecological characteristics of the land surface. From the MODIS archive, we used the average monthly Normalized Difference Vegetation Index (NDVI) to infer vegetation density [40]. NDVI is sometimes also referred to as a measure of “greenness” and is high over forested areas and low over areas with sparse vegetation. In addition, we used the vegetation continuous field [41] product as a measure of the percentage of tree cover. In contrast to NDVI, percent tree cover specifically measures the land surface area covered by trees, disregarding other types of vegetation. From QuickScat (QSCAT), we obtained raw backscatter measurements for the month of August that capture attributes related to surface moisture [40], and, from the Shuttle Radar Topography Mission (SRTM), we acquired elevation data.Time series of the remote sensing data sources were acquired from 2001 for QSCAT and tree cover and means over the period 2000–2004 for NDVI. Variables with native resolutions higher (e.g., SRTM: 30 m) or lower (e.g., QSCAT: 2.25 km) than 1 km were reaggregated to a 1 km grid cell resolution. This resolution is often used in ecological niche modeling at the regional or subcontinental scales and balances resolutions from climate data, which often have coarser native resolutions, and remote sensing data, with higher native resolutions. To improve interpretation, we checked for covariance among variables and only included those with substantial unique variance. Various criteria were used to decide which layers of correlated pairs (with Pearson's correlations on the order of 0.9 or larger) were retained for further analysis. These included keeping layers that are more commonly used in distribution modeling (WorldClim) or that exhibit larger contrast/variance over the study area (QSCAT) as well as having best data quality (NDVI). Pearson’s cross correlations of the used environmental variables are shown in supplementary Table S2.


## 2.2.2. Generalized Dissimilarity Modeling

To analyze the relative contribution of geographic distance and environmental variables to explaining the similarity among sampling locations in the percentage of hybrids, we used generalized dissimilarity modeling (GDM, [42]). GDM is a matrix regression technique that predicts biotic dissimilarity (turnover) between sites based upon environmental heterogeneity and geographic distance. One advantage of GDM over other modeling methodologies is that it can fit nonlinear relationships between predictor and response variables through the use of I-spline basis functions [42]. Another advantage of GDM is that, because it uses pairwise comparisons between sites, it can explicitly take into account the influence of geographic distance on explaining biological variation—a particularly important feature in our study. The relative importance of predictor variables in a GDM can be assessed by means of response curves. To further evaluate the extent to which geographic distance is potentially correlated with environmental differences, we ran independent tests with the following sets of predictor variables: (1) environmental variables and geographic distance; (2) only geographic distance; (3) only environmental variables.GDM is a two-step method: first, dissimilarities of a set of predictor variables are fitted to the genetic dissimilarities (the response variables). The contributions of predictor variables to explaining the observed response variation are tested by permutations, and only those variables that are significant are retained in the final model. These procedures result in a function that describes the relationship between environmental and response variables. Second, using the function resulting from the first step, a spatial prediction is made of the response variable patterns. For visualization purposes, classes of similar response are color coded, where larger color differences between two localities represent larger differences in the proportion of hybrid individuals. In order to assess the significance of the level of variation that was explained by our models, we ran additional models in which the environmental layers were substituted by layers with random values for each grid cell. The resulting percentage of variation explained was compared to that of the full model. We considered the performance of the full model not significant if it explained an equal amount or less of the total variation than a model with random environmental variables.


## 3. Results

### 3.1. Genetic Analyses

Using multilocus genotypes, the proportion of hybrid individuals in the 21 sympatric and parapatric populations ranged from 0.00 to 0.57 using the “80–20” criterion and 0.00 to 1.00 using the “90–10” criterion (Table1). Across all sympatric and parapatric populations, 20 individuals (“80–20” criterion) and 33 individuals (“90–10” criterion) of the 163 sampled (12.2% and 20.2%, resp.) were classified as hybrids. Focusing on the 20 hybrid individuals classified using the “80–20” criterion, hybrids were found as far west as the Roosevelt National Forest sampling locality in Colorado and as far east as the Carpenter Game Production Area along the Missouri River in South Dakota (Figure 1, Table 1). Using the more liberal hybrid classification criterion (“90–10”), the presence of hybrids extended east to the Newton Hills Game Production Area in South Dakota; the westernmost locality where hybrids were present was unchanged (Figure 1, Table 1).


### 3.2. Spatial Analyses

Generalized dissimilarity models (Figure2) performed well as measured by the total variance explained and performed significantly better than models with random environmental variables. The full models explained 58.9% (“90–20” criterion) and 63.9% (“80–20” criterion) of the total observed variation in the similarity of the level of hybridization (Table 4). Models using only environmental variables explained about 6% less of the observed variation, and models using only geographic distance explained 38.8% (“90–10” criterion) and 39.9% (“80–20” criterion) of the variation. The GDM results based on the absolute difference in mean Q-scores were very similar (not shown) and are not discussed further.Table 4
Generalized dissimilarity modeling results. Shown are the percentages of the total variation observed and important predictor variables for “80–20” and “90–10” criteria for full models (using both geographic distance and environmental variables), models with environmental variables only, models with geographic distance only, and random models.
Model% variance explainedVariables selected*80–20 full63.9dist, Bio12, Bio10, elevation, Bio1880–20 environment only57.4Bio16, elevation, Bio12, Bio1080–20 distance only39.9dist80–20 random env0.690–10 full58.9dist, Bio18, Bio10, elevation, Bio1290–10 environment only53.4elevation, Bio18, Bio10, Bio1290–10 distance only38.8Dist90–10 random env0.4*Variables are shown in order of decreasing importance (Figure3); dist: geographic distance; for an explanation of the bioclimatic variables, see Table 3.Figure 2
Generalized dissimilarity modeling results of the proportion of pureP. amoena individuals, plus one half of the hybrid individuals, in a population. Left two panels: “80–20” criterion; right two panels: “90–10” criterion. Top two panels: full model (including both geographic distance and environmental variables); bottom two panels: models for only environmental variables. Colors indicate the proportion of pure P. amoena individuals, plus one half of the hybrid individuals, in the population (brown/tan 100%, blue = 0%).The most important predictor variable in the full models was geographic distance, whereas these were Bio18 (precipitation of the warmest quarter) and elevation for models using only environmental variables (Figure3). These results suggest that geographic distance and environmental heterogeneity along our sampling sites are partially correlated, but that geographic distance alone does not perform as well in explaining differences in the level of hybridization as environmental variables or the combination of distance and environment.Figure 3
Response curves for generalized dissimilarity models of the proportion of pureP. amoena individuals, plus one half of the hybrid individuals, in a population. Results are shown for “80–20” and “90–10” criteria, and for full models (with both geographic distance and environmental variables) and models using only environmental variables. Response curves in GDM are nonlinear, but constrained to be monotonic. The slope of each function is indicative of the rate of change in the proportion of hybrids along the environmental gradient concerned. The maximum height of a response curve indicates the relative importance of that variable in explaining the observed variation in the  similarity between populations in the proportion of hybrids. The y-axes indicate fitted functions f(x), where x denotes the environmental variable concerned. Each of the functions is fitted as a linear combination of I-spline basis functions [41].


## 3.1. Genetic Analyses

Using multilocus genotypes, the proportion of hybrid individuals in the 21 sympatric and parapatric populations ranged from 0.00 to 0.57 using the “80–20” criterion and 0.00 to 1.00 using the “90–10” criterion (Table1). Across all sympatric and parapatric populations, 20 individuals (“80–20” criterion) and 33 individuals (“90–10” criterion) of the 163 sampled (12.2% and 20.2%, resp.) were classified as hybrids. Focusing on the 20 hybrid individuals classified using the “80–20” criterion, hybrids were found as far west as the Roosevelt National Forest sampling locality in Colorado and as far east as the Carpenter Game Production Area along the Missouri River in South Dakota (Figure 1, Table 1). Using the more liberal hybrid classification criterion (“90–10”), the presence of hybrids extended east to the Newton Hills Game Production Area in South Dakota; the westernmost locality where hybrids were present was unchanged (Figure 1, Table 1).


## 3.2. Spatial Analyses

Generalized dissimilarity models (Figure2) performed well as measured by the total variance explained and performed significantly better than models with random environmental variables. The full models explained 58.9% (“90–20” criterion) and 63.9% (“80–20” criterion) of the total observed variation in the similarity of the level of hybridization (Table 4). Models using only environmental variables explained about 6% less of the observed variation, and models using only geographic distance explained 38.8% (“90–10” criterion) and 39.9% (“80–20” criterion) of the variation. The GDM results based on the absolute difference in mean Q-scores were very similar (not shown) and are not discussed further.Table 4
Generalized dissimilarity modeling results. Shown are the percentages of the total variation observed and important predictor variables for “80–20” and “90–10” criteria for full models (using both geographic distance and environmental variables), models with environmental variables only, models with geographic distance only, and random models.
Model% variance explainedVariables selected*80–20 full63.9dist, Bio12, Bio10, elevation, Bio1880–20 environment only57.4Bio16, elevation, Bio12, Bio1080–20 distance only39.9dist80–20 random env0.690–10 full58.9dist, Bio18, Bio10, elevation, Bio1290–10 environment only53.4elevation, Bio18, Bio10, Bio1290–10 distance only38.8Dist90–10 random env0.4*Variables are shown in order of decreasing importance (Figure3); dist: geographic distance; for an explanation of the bioclimatic variables, see Table 3.Figure 2
Generalized dissimilarity modeling results of the proportion of pureP. amoena individuals, plus one half of the hybrid individuals, in a population. Left two panels: “80–20” criterion; right two panels: “90–10” criterion. Top two panels: full model (including both geographic distance and environmental variables); bottom two panels: models for only environmental variables. Colors indicate the proportion of pure P. amoena individuals, plus one half of the hybrid individuals, in the population (brown/tan 100%, blue = 0%).The most important predictor variable in the full models was geographic distance, whereas these were Bio18 (precipitation of the warmest quarter) and elevation for models using only environmental variables (Figure3). These results suggest that geographic distance and environmental heterogeneity along our sampling sites are partially correlated, but that geographic distance alone does not perform as well in explaining differences in the level of hybridization as environmental variables or the combination of distance and environment.Figure 3
Response curves for generalized dissimilarity models of the proportion of pureP. amoena individuals, plus one half of the hybrid individuals, in a population. Results are shown for “80–20” and “90–10” criteria, and for full models (with both geographic distance and environmental variables) and models using only environmental variables. Response curves in GDM are nonlinear, but constrained to be monotonic. The slope of each function is indicative of the rate of change in the proportion of hybrids along the environmental gradient concerned. The maximum height of a response curve indicates the relative importance of that variable in explaining the observed variation in the  similarity between populations in the proportion of hybrids. The y-axes indicate fitted functions f(x), where x denotes the environmental variable concerned. Each of the functions is fitted as a linear combination of I-spline basis functions [41].


## 4. Discussion

Understanding the selective forces maintaining hybrid zones is an important component of understanding the evolutionary process generating and maintaining biodiversity. If natural selection in response to different environments during periods of allopatry contributed to divergence and speciation between closely related taxa, we might expect patterns of hybridization upon secondary contact to reflect the patterns of local environmental heterogeneity. In this study, we have shown that patterns of hybridization betweenPasserina amoena and Passerina cyaneaare best explained by a model that includes both patterns of environmental variation across the Rocky Mountains/Great Plains ecotone and geographic distance between sampling localities (Figures 2 and 3, Table 4). Although the distance between sampling localities was clearly an important determinant of the differences in hybridization among populations, environmental variation alone was a more powerful predictor. This result held for both criteria we used to assess whether an individual was a hybrid (“80–20” and “90–10”).This pattern supports the findings of Swenson [29], who used a much different approach to explore the relationship between environmental heterogeneity and hybridization in Passerina buntings. Whereas our analyses focused on the patterns of hybridization within and among populations, the unit of hybridization in Swenson’s analysis was the geographic location of an individual determined to be a hybrid. The proportion of hybrid individuals in a population, which was central to our analysis, did not factor into Swenson’s work. Instead, Swenson solely focused on whether or not a hybrid had been recorded in a particular location and did not explore the similarities in the proportion of hybrids among populations as we did. Although the relevant environmental variables identified in our study (BioClim variables Bio10, Bio12, and Bio18, Table 4) were slightly different from the environmental variable with the greatest explanatory power in Swenson’s study (mean annual temperature, BioClim variable Bio1), the variables are correlated (Supplemental Table S2). For example, the correlation between Bio18 (precipitation during the warmest quarter), which had the greatest explanatory power in both of our “environment only” models (Table 4), and mean annual temperature (Bio1) was 0.506 (Supplementary Table S1). Despite the difference in methodologies between the studies, both found that environmental variation predicts patterns of hybridization, which strongly indicates that local environments influence the structure of the Passerina bunting hybrid zone.At present, our data do not allow us to determine what is driving the relationship between hybridization and environmental variation, but we consider three general, not mutually exclusive possibilities. The first is that hybrid individuals may either prefer or avoid particular habitats, the second is that there may be differences in habitat-specific performance, and the third is that mate choice patterns could differ across the hybrid zone.That the fitness of hybrid individuals may be related to environmental conditions is not new as it is central to the environmental gradient selection hypothesis of hybrid zone structure [43] and its close relative, the bounded-superiority hypothesis [14]. In both these hypotheses, the selection pressures faced by hybrid individuals depend on ecological context. It is important to note that hybrid individuals may or may not actively prefer particular habitats. It is possible that individuals might disperse some given distance and that any subsequent fitness consequences in a particular habitat may not be driven by choice.Our results are consistent with the possibility that the environmental variation present across thePasserina hybrid zone somehow influences hybrid fitness. It has been argued that both the oriole and flicker hybrid zones (Icterus bullockii × Icterus galbula and Colaptes auratus auratus × Colaptes auratus cafer, resp.), which are also located in the eastern Rocky Mountains/western Great Plains ecotone, are maintained because hybrid individuals are restricted to particular habitats. In part due to differences in metabolic performance between the two species, Rising [44] suggested that I. galbula is better adapted to the more mesic environment characteristic of eastern deciduous forest and I. bullockii is better adapted to the more xeric environments that are more dominant in western North America. He further suggested that hybrids could only survive in the intervening ecotone [44]. Similarly, hybridization between the red-shafted and yellow-shafted subspecies of the Northern Flicker (C. auratus) appears to be strongly related to moisture patterns present across the same ecotone [14, 45]. Since mean temperature and precipitation of the warmest quarter (Bio10, Bio18) were important contributors in the GDM results (Table 4), it is possible that P. cyanea and P. amoena are differentially adapted to their respective habitats (more mesic for P. cyaneaxeric for P. amoena) similar to what Rising hypothesized for Icterus [44]. Either through metabolic differences or some other mechanism, hybrid Passerina buntings may suffer greater fitness consequences on one or both sides of the hybrid zone. Irrespective of the exact mechanism, environmentally influenced differences in hybrid fitness suggest that natural-selection-driven adaptations to different habitats during the course of divergence between P. amoenaand P. cyaneacontributed to the speciation process, a hallmark of ecological speciation.The third possible driver of the differences in hybridization relates to mate choice, which may depend on ecological context. If sensory bias plays a role in female mate choice [46] and the local environment influences a female ability to distinguish between con and heterospecific mates, the rate of hybridization may differ among populations. Although natal dispersal patterns are unknown for Passerina buntings [47, 48], if they do not differ among populations, then the relationship between environmental heterogeneity and hybridization observed could potentially result from differences in female mate choice. Previous work indicated that female buntings can readily, but not perfectly, distinguish between con and heterospecific males [19, 20], but it is not known whether the local environment can influence female mate choice abilities.Recent work has shown a significant narrowing and westward shift in the structure of thePasserina hybrid zone over the past 40–45 years [30]. Although climate data from WorldClim [38] are not of sufficient resolution to test this hypothesis, it is intriguing to speculate that changes in precipitation over the same time frame played an important role in the hybrid zone movement. Additionally, climate models suggest that annual precipitation will increase over much of North America, with the exception of the southwest, over the next 100 years [49]. Interestingly, the majority of the increase in annual precipitation is predicted to be driven by precipitation in the winter months (December, January, and February), whereas precipitation in the summer months (June, July, and August) is predicted to decrease [49]. If the structure of the Passerina hybrid zone is mediated through precipitation patterns during the warmest part of the year, as our data suggests, the zone may actually reverse course and shift eastward in the future as the ecotone becomes drier.Regardless of the exact mechanism, it is clear that in order to fully understand hybridization and continued reproductive isolation betweenP. amoena and P. cyaneathe ecological context within which the interactions between the species occur should not be ignored. Given that these species did not diverge sympatrically [22], they are likely adapted to different environmental conditions and these differences influence hybridization patterns. Consequently, the long-term stability of the hybrid zone, as well as the species themselves, depends at least in part on the stability of the patterns of environmental heterogeneity found in the eastern Rocky Mountains/western Great Plains ecotone.



---

*Source: 295463-2012-02-09.xml*