International Journal of EcologyVolume 2012 (2012), Article ID 295463, 11 pagesdoi:10.1155/2012/295463

Research Article

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

Matthew D. Carling¹ and Henri A. Thomassen²

¹Berry Biodiversity Conservation Center, Department of Zoology and Physiology, University of Wyoming, 1000 E. University Avenue, Laramie, WY 82071, USA²Institute of Evolution and Ecology, University of Tübingen, Auf der Morgenstelle 28, 72076 Tübingen, Germany

Received 25 June 2011; Revised 8 September 2011; Accepted 27 September 2011

Academic Editor: Zachariah Gompert

Copyright © 2012 Matthew D. Carling and Henri A. Thomassen. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 between Passerina 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.

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)

(b)

Figure 1: (a) Primary breeding distributions and allopatric sampling localities of Passerina 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.

More recently, coalescent analyses of the evolutionary history of P. 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 of P. 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.

Code
Locality
Latitude
Longitude

Proportion of hybrids (80–20)
Proportion of hybrids (90–10)

1
WA
47.8
−118.9
29
—
—

2
ID
44.1
−116.2
3
—
—

3
MT: Custer National Forest no. 1
45.1
−108.5
4
0.00
0.00

4
WY: Bighorn National Forest no. 1
44.8
−107.3
8
0.00
0.00

5
WY: Bighorn National Forest no. 2
44.6
−107.7
2
0.00
0.00

6
CO: Roosevelt National Forest
40.7
−105.2
7
0.14
0.14

7
MT: Custer National Forest no. 2
45.7
−106
5
0.00
0.00

8
WY: Medicine Bow National Forest
42.4
−105.3
3
0.00
0.00

9
ND: Little Missouri National Grassland
46.8
−103.5
13
0.00
0.00

10
SD: Custer National Forest
45.3
−103.2
10
0.00
0.00

11
WY: Sand Creek*
44.5
−104.1
14
0.36
0.50

12
NE: White River*
42.6
−103.5
8
0.38
0.50

13
SD: Black Hills National Forest
43.7
−103.8
5
0.20
0.20

14
NE: Ponderosa State Wildlife Management Area
42.6
−103.3
7
0.57
1.00

15
SD: The Nature Conservancy Whitney Preserve
43.3
−103.6
3
0.00
0.00

16
NE: Nebraska National Forest
42.8
−102.9
8
0.25
1.00

17
SD: Ft. Meade National Recreation Area
44.4
−103.5
4
0.50
0.50

18
NE: Nenzel*
42.8
−101.1
4
0.25
0.25

19
NE: The Nature Conservancy
42.8
−100
18
0.00
0.00

    Niobrara Valley Preserve

20
SD: Carpenter Game Production Area
43.7
−99.5
13
0.08
0.08

21
ND: The Nature Conservancy
46.3
−97.3
4
0.00
0.00

    Pigeon Point Preserve

22
NE: Wiseman State Wildlife Area
42.8
−97.1
10
0.00
0.00

23
SD: Newton Hills Game Production Area
43.2
−96.6
13
0.00
0.08

24
MN
44.9
−93.7
10
—
—

25
IL
40.5
−88.9
10
—
—

26
MI
42.3
−83.7
7
—
—

*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. cyanea haplotypes [7, 22]. We then generated multilocus genotypes for each individual by coding each allele at each locus as either P. amoena or 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. cyanea individuals.

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 where PL₁ is the proportion of P. amoena individuals in population 1, PH₁ is the proportion of hybrids in population 1, PL₂ is the proportion of P. amoena individuals in population 2, and PH₂ 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 -scores, as estimated by STRUCTURE [35–38].

Table 2: Matrix of hybridization distances ; 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.1
WY: Bighorn NF no.1
WY: Bighorn NF no.2
CO: Roosevelt NF
MT: Custer NF no.2
WY: Medicine Bow NF
ND: Little Missouri NG
SD: Custer NF
WY: Sand Creek
NE: White River
SD: Black Hills NF
NE: Ponderosa SWMA
SD: TNC Whitney Preserve
NE: Nebraska NF
SD: Ft. Meade NRA
NE: Nenzel
NE: TNC Niobrara Valley Preserve
SD: Carpenter GPA
ND: TNC Pigeon Point Preserve
NE: Wiseman SWA
SD: Newton Hills GPA

MT: Custer NF no. 1
—
0.00
0.00
0.07
0.00
0.00
0.00
0.00
0.32
0.19
0.10
0.29
0.00
0.13
0.50
0.88
1.00
0.96
1.00
1.00
1.00

WY: Bighorn NF no. 1
0.00
—
0.00
0.07
0.00
0.00
0.00
0.00
0.32
0.19
0.10
0.29
0.00
0.13
0.50
0.88
1.00
0.96
1.00
1.00
1.00

WY: Bighorn NF no. 2
0.00
0.00
—
0.00
0.00
0.00
0.00
0.00
0.32
0.19
0.10
0.29
0.00
0.13
0.50
0.88
1.00
0.96
1.00
1.00
1.00

CO: Roosevelt NF
0.07
0.07
0.07
—
0.07
0.07
0.07
0.07
0.25
0.12
0.03
0.21
0.07
0.05
0.43
0.80
0.93
0.89
0.93
0.93
0.93

MT: Custer NF no. 2
0.00
0.00
0.00
0.07
—
0.00
0.00
0.00
0.32
0.19
0.10
0.29
0.00
0.13
0.50
0.88
1.00
0.96
1.00
1.00
1.00

WY: Medicine Bow NF
0.00
0.00
0.00
0.07
0.00
—
0.00
0.00
0.32
0.19
0.10
0.29
0.00
0.13
0.50
0.88
1.00
0.96
1.00
1.00
1.00

ND: Little Missouri NG
0.00
0.00
0.00
0.07
0.00
0.00
—
0.00
0.32
0.19
0.10
0.29
0.00
0.13
0.50
0.88
1.00
0.96
1.00
1.00
1.00

SD: Custer NF
0.00
0.00
0.00
0.07
0.00
0.00
0.00
—
0.32
0.19
0.10
0.29
0.00
0.13
0.50
0.88
1.00
0.96
1.00
1.00
1.00

WY: Sand Creek
0.32
0.32
0.32
0.25
0.32
0.32
0.32
0.32
—
0.13
0.22
0.04
0.32
0.20
0.18
0.55
0.68
0.64
0.68
0.68
0.68

NE: White River
0.25
0.25
0.25
0.18
0.25
0.25
0.25
0.25
0.07
—
0.09
0.10
0.19
0.06
0.31
0.69
0.81
0.77
0.81
0.81
0.81

SD: Black Hills NF
0.10
0.10
0.10
0.03
0.10
0.10
0.10
0.10
0.22
0.15
—
0.19
0.10
0.03
0.40
0.78
0.90
0.86
0.90
0.90
0.90

NE: Ponderosa SWMA
0.50
0.50
0.50
0.43
0.50
0.50
0.50
0.50
0.18
0.25
0.40
—
0.29
0.16
0.21
0.59
0.71
0.68
0.71
0.71
0.71

SD: TNC Whitney Preserve
0.00
0.00
0.00
0.07
0.00
0.00
0.00
0.00
0.32
0.25
0.10
0.50
—
0.13
0.50
0.88
1.00
0.96
1.00
1.00
1.00

NE: Nebraska NF
0.50
0.50
0.50
0.43
0.50
0.50
0.50
0.50
0.18
0.25
0.40
0.00
0.50
—
0.38
0.75
0.88
0.84
0.88
0.88
0.88

SD: Ft. Meade NRA
0.50
0.50
0.50
0.43
0.50
0.50
0.50
0.50
0.18
0.25
0.40
0.00
0.50
0.00
—
0.38
0.50
0.46
0.50
0.50
0.50

NE: Nenzel
0.88
0.88
0.88
0.80
0.88
0.88
0.88
0.88
0.55
0.63
0.78
0.38
0.88
0.38
0.38
—
0.13
0.09
0.13
0.13
0.13

NE: TNC Niobrara Valley Preserve
1.00
1.00
1.00
0.93
1.00
1.00
1.00
1.00
0.68
0.75
0.90
0.50
1.00
0.50
0.50
0.13
—
0.04
0.00
0.00
0.00

SD: Carpenter GPA
0.96
0.96
0.96
0.89
0.96
0.96
0.96
0.96
0.64
0.71
0.86
0.46
0.96
0.46
0.46
0.09
0.04
—
0.04
0.04
0.04

ND: TNC Pigeon Point Preserve
1.00
1.00
1.00
0.93
1.00
1.00
1.00
1.00
0.68
0.75
0.90
0.50
1.00
0.50
0.50
0.13
0.00
0.04
—
0.00
0.00

NE: Wiseman SWA
1.00
1.00
1.00
0.93
1.00
1.00
1.00
1.00
0.68
0.75
0.90
0.50
1.00
0.50
0.50
0.13
0.00
0.04
0.00
—
0.00

SD: Newton Hills GPA
0.96
0.96
0.96
0.89
0.96
0.96
0.96
0.96
0.64
0.71
0.86
0.46
0.96
0.46
0.46
0.09
0.04
0.00
0.04
0.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 (Table 3). 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 Record
Instrument
Ecological attributes
Variables derived

Normalized Difference Vegetation Index (NDVI)
Satellite-MODIS
Vegetation density
NDVI mean

Vegetation Continuous Field (VCF)*§
Satellite-MODIS
Forest cover + heterogeneity
Tree cover

Scatterometer-Backscatter*†
Satellite-QSCAT
Surface moisture + roughness (forest structure), seasonality
QSCAT aug

DEM*
SpaceShuttle SRTM
Topography + ruggedness
SRTM

WorldClim¶
Station-Network
Bioclimatic variables
Bio1, 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 (Table 1). 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 (Figure 2) 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 explained
Variables selected*

80–20 full
63.9
dist, Bio12, Bio10, elevation, Bio18

80–20 environment only
57.4
Bio16, elevation, Bio12, Bio10

80–20 distance only
39.9
dist

80–20 random env
0.6

90–10 full
58.9
dist, Bio18, Bio10, elevation, Bio12

90–10 environment only
53.4
elevation, Bio18, Bio10, Bio12

90–10 distance only
38.8
Dist

90–10 random env
0.4

*Variables are shown in order of decreasing importance (Figure 3); dist: geographic distance; for an explanation of the bioclimatic variables, see Table 3.

Figure 2: Generalized dissimilarity modeling results of the proportion of pure P. 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 (Figure 3). 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 pure P. 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 -axes indicate fitted functions , where 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 between Passerina amoena and Passerina cyanea are 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 the Passerina 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. cyanea xeric 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. amoena and P. cyanea contributed 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 the Passerina 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 between P. amoena and P. cyanea the 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.

Acknowledgments

The authors thank S. Birks (University of Washington Burke Museum of Natural History and Culture), A. Capparella (Illinois State University), J. Hinshaw (University of Michigan Museum of Zoology), and M. Westberg (University of Minnesota Bell Museum) for generously providing tissue samples for this project. Scientific collecting permits were granted by the Colorado Division of Wildlife, Montana Fish, Wildlife and Parks, Nebraska Game and Parks, North Dakota Game and Fish, South Dakota Game, Fish and Parks, Wyoming Game and Fish, and the United States Fish and Wildlife Service. Louisiana State University Institutional Animal Care and Use Committee Protocol: 08-025. The map data were provided by NatureServe in collaboration with Robert Ridgely, James Zook, The Nature Conservancy Migratory Bird Program, Conservation International-CABS, World Wildlife Fund, USA, and Environment Canada-WILDSPACE. This work was supported by grants from NSF (DEB-0543562 and DBI-0400797 to R. T. Brumfield at Louisiana State University, DEB-0515981 and DEB-0814277 to I. J. Lovette at Cornell University, and DEB-0808464 to M. D. Carling), as well as by grants from the AMNH Chapman Fund, AOU, Explorer’s Club, Louisiana State University Museum of Natural Science Birdathon and Prepathon Funds, Louisiana State University Department of Biological Sciences, and Sigma-Xi to M. D. Carling.

References

G. M. Hewitt, “Hybrid zones-natural laboratories for evolutionary studies,” Trends in Ecology and Evolution, vol. 3, no. 7, pp. 158–167, 1988.

R. G. Harrison, “Hybrid zones: windows on the evolutionary process,” in Oxford Surveys in Evolutionary Biology, D. J. Futuyma and J. Antonovics, Eds., pp. 68–129, Oxford University Press, Oxford, UK, 1990.

R. G. Harrison, “Hybrids and hybrid zones: historical perspective,” in Hybrid Zones and the Evolutionary Process, R. G. Harrison, Ed., pp. 3–12, Oxford University Press, New York, NY, USA, 1993.

M. L. Arnold, Natural Hybridization and Evolution, Oxford University Press, New York, NY, USA, 1997.

J. Mallet, “Hybridization as an invasion of the genome,” Trends in Ecology and Evolution, vol. 20, no. 5, pp. 229–237, 2005.

B. A. Payseur, J. G. Krenz, and M. W. Nachman, “Differential patterns of introgression across the X chromosome in a hybrid zone between two species of house mice,” Evolution, vol. 58, no. 9, pp. 2064–2078, 2004.

M. D. Carling and R. T. Brumfield, “Speciation in Passerina buntings: introgression patterns of sex-linked loci identify a candidate gene region for reproductive isolation,” Molecular Ecology, vol. 18, no. 5, pp. 834–847, 2009.

O. Seehausen, G. Takimoto, D. Roy, and J. Jokela, “Speciation reversal and biodiversity dynamics with hybridization in changing environments,” Molecular Ecology, vol. 17, no. 1, pp. 30–44, 2008.

M. E. Maan and O. Seehausen, “Ecology, sexual selection and speciation,” Ecology Letters, vol. 14, no. 6, pp. 591–602, 2011.

D. Schluter, “Ecology and the origin of species,” Trends in Ecology and Evolution, vol. 16, no. 7, pp. 372–380, 2001.

J. A. Coyne and H. A. Orr, Speciation, Sinauer Associates, Sunderland, Mass, USA, 2004.

P. Nosil, L. J. Harmon, and O. Seehausen, “Ecological explanations for (incomplete) speciation,” Trends in Ecology and Evolution, vol. 24, no. 3, pp. 145–156, 2009.

J. A. Ender, Geographic Variation, Speciation and Clines, Princeton University Press, Princeton, NJ, USA, 1977.

W. S. Moore, “Evaluation of narrow hybrid zones in vertebrates,” Quarterly Review of Biology, vol. 52, pp. 263–277, 1977.

M. D. Carling and R. T. Brumfield, “Integrating phylogenetic and population genetic analyses of multiple loci to test species divergence hypotheses in Passerina buntings,” Genetics, vol. 178, no. 1, pp. 363–377, 2008.

C. G. Sibley and L. L. J. Short, “Hybridization in the buntings (Passerina) of the Great Plains,” Auk, vol. 76, pp. 443–463, 1959.

S. T. Emlen, J. D. Rising, and W. L. Thompson, “Behavioral and morphological study of sympatry in Indigo and Lazuli buntings of the Great Plains,” Wilson Bulletin, vol. 87, pp. 145–179, 1975.

R. L. Kroodsma, “Hybridization in buntings (Passerina) in North Dakota and eastern Montana,” Auk, vol. 92, pp. 66–80, 1975.

M. C. Baker and A. E. M. Baker, “Reproductive behavior of female buntings: isolating mechanisms in a hybridizing pair of species,” Evolution, vol. 44, no. 2, pp. 332–338, 1990.

M. C. Baker, “Female buntings from hybridizing populations prefer conspecific males,” Wilson Bulletin, vol. 108, no. 4, pp. 771–775, 1996.

M. C. Baker and J. T. Boylan, “Singing behavior, mating associations and reproductive success in a population of hybridizing Lazuli and Indigo Buntings,” Condor, vol. 101, no. 3, pp. 493–504, 1999.

M. D. Carling, I. J. Lovette, and R. T. Brumfield, “Historical divergence and gene flow: coalescent analyses of mitochondrial, autosomal and sex-linked loci in Passerina buntings,” Evolution, vol. 64, no. 6, pp. 1762–1772, 2010.

M. D. Carling and R. T. Brumfield, “Haldane's rule in an avian system: using cline theory and divergence population genetics to test for differential introgression of mitochondrial, autosomal, and sex-linked loci across the Passerina bunting hybrid zone,” Evolution, vol. 62, no. 10, pp. 2600–2615, 2008.

J. B. S. Haldane, “Sex ratio and unisexual sterility in hybrid animals,” Journal of Genetics, vol. 12, no. 2, pp. 101–109, 1922.

O. M. Ocon-Grove, S. Maddineni, G. L. Hendricks, et al., “Pituitary progesterone receptor expression and plasma gonado- tropic concentrations in the reproductively dysfunctional mutant restricted ovulator chicken,” Domestic Animal Endocrinology, vol. 32, pp. 201–210, 2007.

T. Dobzhansky, Genetics and the Origin of Species, Columbia University Press, New York, NY, USA, 1937.

H. J. Muller, “Bearing of the Drosophila work on systematics,” in The New Systematics, J. B. S. Haldane, Ed., pp. 185–268, Clarendon Press, Oxford, UK, 1940.

H. J. Muller, “Isolating mechanisms, evolution, and temperature,” Biological Symposia, vol. 6, pp. 71–125, 1942.

N. G. Swenson, “Gis-based niche models reveal unifying climatic mechanisms that maintain the location of avian hybrid zones in a North American suture zone,” Journal of Evolutionary Biology, vol. 19, no. 3, pp. 717–725, 2006.

M. D. Carling and B. Zuckerberg, “Spatio-temporal changes in the genetic structure of the Passerina bunting hybrid zone,” Molecular Ecology, vol. 20, no. 6, pp. 1166–1175, 2011.

R. S. Ridgely, T. F. Allnut, T. Brooks, et al., Digital Distribution Maps of the Birds of the Western Hemisphere, Version 1.0, NatureServe, Arlington, Va, USA, 2003.

NatureServe, “NatureServe explorer: an online encyclopedia of life, version 6.1,” 2010, http://www.natureserve.org/explorer

M. Stephens, N. J. Smith, and P. Donnelly, “A new statistical method for haplotype reconstruction from population data,” American Journal of Human Genetics, vol. 68, no. 4, pp. 978–989, 2001.

M. Stephens and P. Donnelly, “A comparison of Bayesian methods for haplotype reconstruction from population genetic data,” American Journal of Human Genetics, vol. 73, no. 5, pp. 1162–1169, 2003.

J. K. Pritchard, M. Stephens, and P. Donnelly, “Inference of population structure using multilocus genotype data,” Genetics, vol. 155, no. 2, pp. 945–959, 2000.

D. Falush, M. Stephens, and J. K. Pritchard, “Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies,” Genetics, vol. 164, no. 4, pp. 1567–1587, 2003.

D. Falush, M. Stephens, and J. K. Pritchard, “Inference of population structure using multilocus genotype data: dominant markers and null alleles,” Molecular Ecology Notes, vol. 7, no. 4, pp. 574–578, 2007.

M. J. Hubisz, D. Falush, M. Stephens, and J. K. Pritchard, “Inferring weak population structure with the assistance of sample group information,” Molecular Ecology Resources, vol. 9, no. 5, pp. 1322–1332, 2009.

R. J. Hijmans, S. E. Cameron, J. L. Parra, P. G. Jones, and A. Jarvis, “Very high resolution interpolated climate surfaces for global land areas,” International Journal of Climatology, vol. 25, no. 15, pp. 1965–1978, 2005.

R. B. Myneni, S. Hoffman, Y. Knyazikhin, et al., “Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data,” Remote Sensing of Environment, vol. 83, no. 1-2, pp. 214–231, 2002.

M. C. Hansen, R. S. DeFries, J. R. G. Townshend, R. Sohlberg, C. Dimiceli, and M. Carroll, “Towards an operational MODIS continuous field of percent tree cover algorithm: examples using AVHRR and MODIS data,” Remote Sensing of Environment, vol. 83, no. 1-2, pp. 303–319, 2002.

S. Ferrier, G. Manion, J. Elith, and K. Richardson, “Using generalized dissimilarity modelling to analyse and predict patterns of beta diversity in regional biodiversity assessment,” Diversity and Distributions, vol. 13, no. 3, pp. 252–264, 2007.

R. M. May, J. A. Endler, and R. E. McMurtrie, “Gene frequency clines in the presence of selection opposed by gene flow,” American Naturalist, vol. 109, pp. 659–676, 1975.

J. D. Rising, “A comparison of metabolism and evaporative water loss of baltimore and bullock orioles,” Comparative Biochemistry And Physiology, vol. 31, no. 6, pp. 915–925, 1969.

W. S. Moore and J. T. Price, “Nature of selection in the Northern Flicker hybrid zone and its implications for speciation theory,” in Hybrid Zones and the Evolutionary Process, R. G. Harrison, Ed., pp. 196–225, Oxford University Press, New York, NY, USA, 1993.

M. J. Ryan and A. Keddy-Hector, “Directional patterns of female mate choice and the role of sensory biases,” American Naturalist, vol. 139, pp. S4–S35, 1992.

E. Greene, V. R. Muehter, and W. Davison, “Lazuli Bunting (Passerina amoena),” in The Birds of North America, A. Poole, Ed., Cornell Lab of Ornithology, Ithaca, NY, USA, 1996.

R. B. Payne, “Lazuli Bunting (Passerina amoena),” in The Birds of North America, A. Payne, Ed., 2006.

S. Solomon, D. Qin, M. Manning, et al., Eds., Contribution of Working Group 1 to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, New York, NY, USA, 2007.