Species range expansions facilitated by global climate change have been documented across many taxa. The ecological and evolutionary costs of range expansion in response to climate change are beginning to be teased apart, and have the potential to be strikingly different among taxa experiencing different types of range expansion across highly variable landscapes. We investigated how population expansion and connectivity have affected genetic diversity in the Virginia's Warbler (Oreothlypis virginiae), a recent colonist of the Black Hills of South Dakota, USA. To investigate population connectivity, we sampled Virginia's Warbler across their breeding range. Genetic data from the mtDNA NADH dehydrogenase subunit 2 (ND2) gene and 7 polymorphic microsatellite loci were used to characterize structure within and among populations and to investigate genetic variability associated with connectivity. The analyses suggested little genetic differentiation, minimal population structure, and similar mtDNA haplotype and nuclear heterozygosity diversities throughout all 7 of the sampled regions. Results from Tajima's D and Fu's FS neutrality tests and a starburst haplotype network indicated a demographic expansion similar to what would be expected in a species that underwent historical range expansion following the last glacial maximum. This study demonstrates that range expansion with recurring gene flow can curtail the loss of genetic diversity and prevent significant differentiation in newly colonized areas.
INTRODUCTION
Climate change has naturally occurred throughout Earth's history. However, the current interglacial period, the Holocene, has experienced unexpected warming patterns over the last few thousand years (IPCC 2014). These unexpected patterns are characterized by the rapid speed at which warming trends are taking place. From 1880 to 2012, the global temperature increased by 0.85°C (IPCC 2014). Consequently, global climate change is beginning to have a direct impact on the distributions of contemporary species (Chen et al. 2011), especially in areas of profound climatic shifts (Jansson and Dynesius 2002). Patterns of biodiversity loss may become apparent due to the merging of gene pools and the inability of specialists to adapt to changing environmental conditions (Saccheri et al. 1998, Scriber 2011). Populations incapable of adjusting or species with restricted dispersal capabilities may become isolated in refugia, and in some cases may become extinct (Holt 1990, Saccheri et al. 1998, Parmesan 2006).
Some species adjust to changing conditions by tracking their climatic niche (Walther et al. 2002, Parmesan 2006). Species' ranges are highly dynamic, and can expand, contract, and otherwise change over time (Sexton et al. 2009). Range limits are determined by biotic factors, such as intraspecific or interspecific interactions, and abiotic factors, including climate (Parmesan 2006, Angert 2009). Variations in these biotic and abiotic factors determine the degree of gene flow between and among populations and ultimately influence biodiversity (Slatkin 1987, Pimm 2008). Accordingly, the ecological and evolutionary costs of range expansion in response to climate change are beginning to be teased apart, and have the potential to be strikingly different among taxa experiencing different types of range expansion across highly variable landscapes. Most species' ranges show patterns of expansion at the leading edge associated with shifts in latitude and elevation (Walther et al. 2002, Parmesan 2006). For instance, distributional shifts facilitated by global climate change appear to be moving toward the poles at an average rate of ~6.1 km per decade (Parmesan and Yohe 2003).
Range expansion with recurring gene flow could curtail the loss of neutral and functional genetic diversity in newly colonized areas. In this case, genetic diversity is maintained in peripheral populations through immigration and multiple introduction events (Bialozyt et al. 2005). Alternatively, range expansion across uninhabitable regions to habitat islands could produce founder events, and without recurring gene flow could substantially decrease genetic variation in newly colonized areas. Over time, a sustained lack of gene flow could lead to pronounced genetic structure among populations within a species (Bohonak 1999). When a range shift occurs across a habitat gap, the viability of the newly established population is likely related to the frequency of migration of conspecifics into the new area (Hansson et al. 2000) and the size of the founding population (Sakai et al. 2001). Founder events and genetic bottlenecks resulting from distributional shifts affect the sustainability and persistence of populations by causing a reduction in genetic variation (Falconer 1981, Slatkin 1987). Moreover, without adequate connectivity, genetic structure can arise quickly at the periphery of the expanding population, where the number of individuals leading the expansion can be quite low (Wade and McCauley 1988, Ibrahim et al. 1996, Austerlitz et al. 1997, Le Corre and Kremer 1998).
While it is important to ascertain the effects of distributional shifts associated with human-induced climate change, it can be difficult to differentiate between these shifts and shifts associated with climate change since the last glacial maximum. During the last glacial maximum, many species became isolated in refugia, but as the glacial ice sheets retreated and species dispersed, genetic differentiation and divergence ensued (Hewitt 2000). Following glacial retreat, genetic diversity within a species is often lost through range expansion and subsequent founder events (Hewitt 1999), with the degree of loss positively correlated with the geographic distance from the refugium. Similar patterns of loss of neutral genetic diversity should be found in contemporary range shifts, but since current climate changes are taking place more rapidly, the genetic changes occurring in populations undergoing current expansions should look aberrant when compared with populations in the center of a species' current distribution. Thus far, relatively few empirical studies have investigated the genetic responses of populations undergoing contemporary range expansions or distributional shifts (Garroway et al. 2011, Song et al. 2013, Swaegers et al. 2014).
The Virginia's Warbler (Oreothlypis virginiae) serves as an ecologically interesting species to study the effects of a contemporary range expansion. The Virginia's Warbler is a migratory and philopatric species, with a summer breeding range concentrated in the southwestern United States that extends from southeastern Idaho into mountainous regions of southern Arizona, New Mexico, and western Texas (Olson and Martin 1999). Disjunct populations are located in the Black Hills of South Dakota, central Wyoming, northern Colorado, and southeastern Idaho (Dunn and Garrett 1997). The Black Hills (South Dakota) population of Virginia's Warbler is of particular interest because it is the result of a recent colonization event that occurred ~18 yr ago at the species' northernmost range limit (first occurrence documented in July of 1997; Palmer 1998). The colonization time has been estimated from annual breeding bird surveys conducted since the 1980s.
Comprising 15,540 km2, the Black Hills are an isolated group of mountains located primarily in western South Dakota and extending into northeastern Wyoming (Froiland 1999). They were formed as a domal uplift that occurred concurrently with the Laramide orogeny at the start of the Cenozoic, with subsequent volcanic activity in the northern Black Hills establishing their final shape 37 mya (Froiland 1999). The forests of the Black Hills are surrounded by the mixed-grass prairie of the western Great Plains. Prairie habitat is not suitable for Virginia's Warbler, which occupy dense, shrubby vegetation such as mountain mahogany (Cercocarpus montanus), with interspersed taller vegetation such as ponderosa pine (Pinus ponderosa) and juniper (Juniperus spp.; Dunn and Garrett 1997). Approximately 100 breeding pairs of Virginia's Warbler are summer residents in the mountain mahogany–dominated vegetation community of the southern Black Hills (Swanson et al. 2000). However, the species' total range and current relative abundance in the Black Hills have not yet been determined and require further study. Currently, the closest breeding population of Virginia's Warbler is ~370 km away in the Laramie Mountains of Wyoming.
Records from annual breeding bird surveys suggest that the Virginia's Warbler is experiencing a range expansion rather than a range shift (Sauer et al. 2014). For example, recent sightings of Virginia's Warblers have been reported in southeastern Texas, part of the species' southernmost breeding range. These observations, provided by the North American Breeding Bird Survey (Sauer et al. 2014), demonstrate that the southernmost population is doing well. Dedicated ornithological citizen science organizations, including eBird (ebird.org), corroborate this observation.
Undoubtedly, species capable of adapting to changing conditions and expanding their ranges will be faced with variations in climate and complex topographical features. These factors will determine the degree of gene flow, and ultimately will define species' distributions. Therefore, the goal of this study was to understand how connectivity shapes genetic structure during a contemporary range expansion across a dynamic landscape by testing the following hypotheses: (1) From patterns observed during annual breeding bird surveys, we predicted that the range expansion exhibited by Virginia's Warbler is ongoing. Therefore, introduction of new alleles from migrants could be supplementing the fragmented populations, eroding genetic structure and any signature of a population bottleneck within the Black Hills population; (2) Given the relatively recent colonization of the Black Hills by Virginia's Warbler, we did not expect to detect any significant genetic differentiation between this population and others. Even if gene flow is low, there have not been enough generations for genetic drift to have had a profound impact on the population; (3) Due to recent colonization and expected gene flow, we did not anticipate a strong signature of isolation-by-distance in Virginia's Warbler populations, but we did expect to observe patterns indicative of a rapid and recent population expansion; and (4) We did not expect to see a profound latitudinal allele frequency gradient along the axis of expansion, as gene flow should erode signatures of founder events during the range expansion process. To test these hypotheses, we analyzed genetic variation (mtDNA and 7 microsatellite loci) in Virginia's Warblers across their range.
METHODS
Sample Collection and DNA Extraction
We captured 20 Virginia's Warblers using mist nets at 3 sites in the Black Hills of South Dakota. All individuals were captured during the breeding season, with fieldwork concentrated in June 2012. Two tail feathers and one secondary feather were collected from each individual. To sample across the breeding range, an additional 112 tissue samples were acquired from several museum collections and included samples from Nevada, New Mexico, Arizona, Colorado, and Utah (Museum of Vertebrate Zoology at Berkeley, n = 5; Museum of Southwest Biology, n = 12; Burke Museum, n = 23; Louisiana State University Museum of Zoology, n = 5; Denver Museum of Nature and Science, n = 14; American Museum of Natural History, n = 9; and Utah Museum of Natural History, n = 44). We collected a total of 132 samples from throughout the distribution of the Virginia's Warbler (Figure 1). From overall spatial clustering of individuals and the geographic proximity of most locations sampled, we grouped individuals into 8 separate populations, as follows: Nevada (n = 23); South Dakota (n = 20); New Mexico (n = 14); Arizona (n = 14); Colorado (n = 15); northern Utah (n = 16); southern Utah (n = 26); and Mexico (n = 4). Utah was split into northern and southern regions in an attempt to equalize the number of individuals from each location.
FIGURE 1.
Sampling site locations for 132 individuals of Virginia's Warbler divided into 8 geographic regions: South Dakota (n = 20); Colorado (n = 15); New Mexico (n = 14); Arizona (n = 14); Nevada (n = 23); northern Utah (n = 16); southern Utah (n = 26); and Mexico (n = 4). The Virginia's Warbler breeding range is outlined in short red dashed lines, its proposed migration route is outlined in long black dashed lines, and its overwintering area in southern Mexico is outlined in short blue dashed lines.
Total genomic DNA was extracted from tissue samples using a Qiagen DNeasy tissue extraction kit (Qiagen, Valencia, California, USA). Fifty-eight of the 132 samples were toepad clips from prepared museum specimens, and thus the DNA extraction procedure included modifications as outlined by Irestedt et al. (2006). To avoid contaminating the historical samples, isolation of DNA from toe pads was performed in a separate clean lab (a lab never before exposed to avian DNA and connected to a separate air handling system). DNA from tail and secondary feathers was isolated following a standard phenol:chloroform extraction protocol (Sambrook et al. 1989).
Mitochondrial DNA Sequencing
A total of 1,041 base pairs (bp) of the mtDNA NADH dehydrogenase subunit 2 (ND2) gene was amplified by PCR using established primers L5215 (Hackett 1996) and H6313 (Johnson and Sorenson 1998). To ensure complete amplification of ND2 using toe pad samples, internal primers were designed from aligned Virginia's Warbler ND2 sequences using Primer Express (Applied Biosystems, Foster City, California, USA) to break up the gene into 4 smaller fragments (Appendix Table 4). We were unable to amplify 9 samples: JP 739, JP 741, JP 742, UMNH 8857, UMNH 12819, UMNH 15378, UMNH 19254, UMNH 21323, and UMNH 22357. PCR was performed in 15 μl reactions using a cycling protocol that included denaturation at 95°C, followed by 35 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 45 s. This was followed by a 10 min extension at 72°C and a 4°C soak. Following PCR, products that produced a clear band following electrophoresis on an agarose gel were purified using ExoSAP-IT (USB, Cambridge, Massachusetts, USA). Sequencing included 10 μl BigDye cycle sequencing reactions (Applied Biosystems) using the following parameters: denaturation at 96°C for 60 s, and 25 cycles of 96°C for 10 s, 50°C for 05 s, and 60°C for 90 s. Reactions were purified using Agencourt CleanSEQ magnetic beads (Agencourt Bioscience, Beverly, Massachusetts, USA). Samples were eluted with 0.1 mM EDTA, followed by sequencing on an ABI 3130 genetic analyzer (Applied Biosystems). Sequences were edited and aligned using Sequencher 5.0 (Gene Codes, Ann Arbor, Michigan, USA).
Microsatellite Characterization
Using primers previously isolated and designed for a closely related species, the Golden-winged Warbler (Vermivora chrysoptera; Stenzler et al. 2004), all Virginia's Warbler samples were initially screened at 11 microsatellite loci (Appendix Table 5). A universal M13 sequence was incorporated at the 5′ end of each forward primer, onto which a fluorescent label (FAM-10, HEX, and NED) was then attached during PCR. Microsatellite amplifications were performed in 15 μl reactions following a protocol outlined by Stenzler et al. (2004), but with a melting temperature optimized at 58°C. PCR cycling conditions consisted of an initial denaturation cycle at 94°C for 5 min; 10 cycles of 94°C for 30 s, 58°C for 30 s, and 72°C for 30 s; 30 cycles of 94°C for 30 s, 48°C for 30 s, and 72°C for 30 s; and a final elongation at 72°C for 15 min. Loci were coloaded for fragment analysis. The 1st coloaded set consisted of FAM-10-labeled loci VeCr02, VeCr03, and VeCro4, and NED-labeled locus VeCr05. The 2nd coloaded set included FAM-10-labeled loci VeCr06, VeCr07, and VeCr10, HEX-labeled loci VeCr08 and VeCr16, and NED-labeled locus VeCr14. Fragments were genotyped on an ABI 3130 genetic analyzer using a ROX 500 bp size standard (Applied Biosystems). Genotyping failed to yield results in 3 loci: VeCr06, VeCr14, and VeCr16. Genotypes and allele calls were scored using Genemapper 4.0 (Applied Biosystems). Individuals with missing data at >2 loci were discarded from genetic analysis (a total of 13 individuals: VAWA 13, VAWA 15, UMNH 19250, UMNH 15378, UMNH 22357, UMNH 6391, UMNH 8857, UMNH 12819, UMNH 15230, UMNH 19251, UMNH 19254, UMNH 21323, and UMNH 21324).
Mitochondrial DNA Data Analyses
To analyze DNA sequence variation, genetic diversity indices were generated using program DNA Sequence Polymorphism 5.0 (DnaSP 5.0; Librado and Rozas 2009). The indices calculated included nucleotide and haplotype diversity (2 measures of DNA polymorphism), haplotype number, and number of private haplotypes. To display the distribution of haplotypes among populations of Virginia's Warbler, sequences from 123 individuals were used to generate a median-joining haplotype network using program Network 4.6.1.1 (fluxus-engineering.com; Bandelt et al. 1999). The 4 individuals sampled in Mexico were collected in April during migration and represent overwintering Virginia's Warblers, not breeding birds. However, sequences from these individuals were included in the haplotype network to ensure that these sequences fell within the network. Their inclusion did not affect mtDNA population genetic estimates, and they were removed from most analyses (except where their inclusion is explicitly noted).
We tested for departures from neutrality using Fu's FS and Tajima's D statistics using DnaSP 5.0 (Librado and Rozas 2009). Fu's FS measures selective neutrality and is sensitive to detecting recent population expansion (Fu 1997.) Negative FS values reflect an excess number of alleles, which is indicative of a recent population expansion and/or genetic hitchhiking (Fu 1997). Positive Fu's FS values suggest a stable population. We also used Tajima's D to test whether selection is acting at a specific locus and whether the overall population of Virginia's Warbler is in stasis (Tajima 1989). Negative values indicate population expansion and/or purifying selection, whereas positive values indicate a decrease in population size and/or balancing selection. Time since expansion (t) was calculated by t = τ/2u (Schenekar and Weiss 2011), where u is the mean substitution rate of ND2 (1.23 × 10−8; Smith et al. 2011). DnaSP 5.0 (Librado and Rozas 2009) was used to estimate tau (τ). We also estimated the timing and extent of population expansion using a Bayesian Skyline analysis as implemented in BEAST (Drummond et al. 2005) by pooling sequences from all 123 individuals. We ran 3 independent runs for 40,000,000 generations, sampling trees every 1,000 generations using an HKY model of sequence evolution. Log Combiner (Drummond et al. 2012) was used to combine the log and tree files from the independent runs, and the first 10% of trees were discarded as burn-in. The trace files from the 3 independent runs were combined to ensure that effective sample size (ESS) values for all estimated parameters exceeded 200.
To evaluate population structure, an analysis of molecular variance (AMOVA) was performed in Arlequin 3.5 (Excoffier and Lischer 2010). Pairwise GST values for mtDNA were estimated using DnaSP. Arlequin was used to calculate mitochondrial-based pairwise FST values. The number of permutations was set at 10,000 at P = 0.05. To investigate the relationship between genetic (FST/(1 − FST)) and geographic distances, we tested for isolation-by-distance (IBD) using Isolation By Distance Web Services 3.23 (IBDWS; Jensen et al. 2005). Pairwise geographic distances were determined by measuring the distance between the center of each population. A Mantel test was performed in IBDWS to detect a matrix correlation between genetic distance and geographic distance. A total of 10,000 randomizations was applied.
Microsatellite Data Analyses
To test for the presence of null alleles and scoring errors due to stuttering or large allele dropout, MICRO-CHECKER 2.2.3 (van Oosterhout et al. 2004) was implemented using Bonferroni correction. Program CONVERT 1.31 (Glaubitz 2004) was used to format the dataset for use in various software packages. Summary statistics, including allele frequencies and expected (He) and observed heterozygosities (Ho), were estimated in GenAlEx 5.01 (Peakall and Smouse 2006, 2012). We used the web version of GENEPOP 4.03 (Raymond and Rousset 1995, Rousset 2008) to test for deviations from Hardy-Weinberg equilibrium and linkage disequilibrium. The Markov chain parameters were raised to 10,000 dememorizations, 1,000 batches, and 10,000 iterations for each assessment. To investigate gene flow, NM estimates were generated using GenAlEx.
Recently established populations may be founded by a small number of individuals. In those populations, loss of heterozygosity and rare alleles are signatures of genetic bottlenecks. Because novel alleles are reduced more rapidly than heterozygosity, under expectations of mutation-drift equilibrium, excess heterozygosity is predicted (Luikart and Cornuet 1998). We tested for signatures of a bottleneck in Virginia's Warbler populations using BOTTLENECK 1.2.02 (Piry et al. 1999). As deemed appropriate for microsatellite use and to account for fewer than 20 loci, we conducted the Wilcoxon signed-rank test with the 2-phase mutation model (TPM; 70% stepwise mutations and variance = 30; Di Rienzo et al. 1994, Piry et al. 1999). The allele frequency distribution test was also performed, with the expectation that effects of a bottleneck would create a shift in the distribution of alleles when compared with the frequencies under mutation-drift equilibrium.
An AMOVA was performed in Arlequin 3.5 (Excoffier and Lischer 2010) to describe genetic variation within and among populations of Virginia's Warbler. Population differentiation was evaluated using FST values, which were also generated in Arlequin 3.5. The number of permutations was set at 10,000 with P = 0.05. The parameter allowing for a level of missing data was raised to 0.30 to account for absent values at various loci. Given our use of microsatellites, we used CoDiDi 1.0 (Correlation between Diversity and Differentiation; Wang 2015) to ensure that FST estimates were not underestimating population differentiation. Among other things, CoDiDi is used to test for a correlation between GST and HS values, allowing for confidence in reporting population differentiation statistics (Wang 2015). We ran CoDiDi with 100,000 permutations. Using pairwise Slatkin linearized FST values (FST/(1 − FST); Slatkin 1995) generated in Arlequin, we used Isolation By Distance Web Services 3.23 (IBDWS; Jensen et al. 2005) to investigate a relationship between genetic and geographic distance with 10,000 randomizations. Again, pairwise geographic distances were determined by measuring the distance between the center of each population. A Mantel test was performed to test for a positive correlation between the genetic and geographic distances.
To ensure that the grouping of individuals (by state) had no influence on population estimates, we alternatively grouped individuals (those with microsatellite data) into 6 populations according to major biogeographic features acting as potential barriers to gene flow in the southwestern United States (Great Basin, n = 24; Wasatch Range, n = 16; eastern Rocky Mountains, n = 14; western Rocky Mountains, n = 17; southern mountain ranges of Arizona and New Mexico, n = 22; and the Black Hills, n = 18). An AMOVA was performed and pairwise FST values were generated using Arlequin 3.5 (Excoffier and Lischer 2010), with the number of permutations set at 10,000 and P = 0.05. The level for missing data was set at 0.30.
Bayesian clustering was done with program Structure (Pritchard et al. 2000) to assign individuals to populations or clusters (K). All analyses were performed inferring lambda and with a burn-in period of 100,000, followed by 100,000 Markov chain Monte Carlo iterations as recommended by Gilbert et al. (2012). The number of groups was set to 1–10, with 20 iterations per K. We initially analyzed individuals without a priori population information, but in an attempt to detect subtle genetic structure, a second analysis was performed with the assumption of correlated allele frequencies among populations and population location priors. Results were evaluated for ΔK (Evanno et al. 2005) and the maximum probability of each K (Pritchard et al. 2000) in Cluster Markov Packager Across K (CLUMPAK; Kopelman et al. 2015).
RESULTS
Genetic Diversity
We analyzed all 1,041 bp of the ND2 gene for 123 samples of Virginia's Warbler (GenBank accession numbers KT963321–KT963443). We found 29 polymorphic sites and 30 haplotypes among the 8 populations (including Mexico). One widespread haplotype was common to all populations. Six additional haplotypes were shared by at least 2 populations. Most, however, were private and singleton haplotypes unique to the population in which they were found. Haplotype diversity (h) was moderate to high for each population, varying from h = 0.47 in the Colorado population to h = 0.84 in the New Mexico population, with an overall diversity of 0.61 (Table 1). Nucleotide diversity (π) was low, with a mean of 0.00103, ranging from π = 0.00058 in the Colorado population to π = 0.00160 in the New Mexico population. Colorado was the least diverse population, with the lowest haplotype diversity and the lowest nucleotide diversity.
TABLE 1.
Mitochondrial DNA genetic diversity indices for 7 sampling localities of Virginia's Warbler in the U.S.
We found linkage disequilibrium in 4 tests, with all 4 being distinct pairs of loci (P < 0.05). Physical linkage of loci could therefore be rejected, with confidence that each locus was independent. MICRO-CHECKER revealed the possibility of null alleles at loci VeCr02, VeCr05, and VeCr07. No evidence for scoring error was detected in loci VeCr02 or VeCr05, while VeCr07 may have been subject to scoring error due to stuttering. Because VeCr02 and VeCr05 were each present in only 1 population, these loci were retained in subsequent genetic analysis. Similarly, an excess of homozygotes in VeCr07 was detected in only 3 populations, 2 of which were in Utah. Downstream analysis was performed with and without VeCr07 in the dataset.
Five populations deviated from Hardy-Weinberg equilibrium: Nevada at VeCr05 (P = 0.02) and VeCr07 (P = 0.001); Arizona at VeCr07 (P = 0.04); Colorado at VeCr08 (P = 0.03); northern Utah at VeCr02 (P = 0.004) and VeCr07 (P = 0.07); and southern Utah at VeCr07 (P = 0.00). The FIS estimate for the Colorado population at locus VeCr08 was low (W&C FIS = −0.68 and R&H global FIS = −0.39), and thus was retained. Since no Hardy-Weinberg equilibrium deviations were detected in more than 1 population for loci VeCr02, VeCr05, and VeCr08, we kept these loci for all future analyses.
Microsatellite polymorphism ranged from 4 (both VeCr03 and VeCr04) to 15 (VeCr10) alleles per locus, with a mean of 7.14 alleles per locus. Arizona had the highest number of alleles (5.43 ± 1.11), while Colorado had the lowest number of alleles (4.29 ± 0.84; Table 2). Observed population heterozygosities were moderate and ranged from 0.44 ± 0.11 in Nevada to 0.61 ± 0.11 in New Mexico, with a grand mean of 0.55 ± 0.04 across all loci and populations.
TABLE 2.
Genetic diversity indices from 7 polymorphic microsatellite loci for 7 Virginia's Warbler sampling localities in the U.S.: Number of effective alleles (Ne), observed heterozygosity (Ho), and expected heterozygosity (He). SE = standard error.
Population Structure
Using mtDNA sequence data, an analysis of molecular variance (AMOVA) revealed more variance within populations (99.54%) than among populations (0.46%). No significant differentiation was detected between populations using mitochondrial-based pairwise GST and FST values (Appendix Table 6). Pairwise FST values ranged from −0.038 to 0.056. A Mantel test revealed no significant relationship between genetic (FST/(1 − FST)) and geographic distances (r = −0.004, P = 0.39) for mtDNA (Appendix Figure 3).
An AMOVA performed on the microsatellite data showed little variance among populations (0.34%). More variance was detected within populations (99.66%). Microsatellite pairwise FST values ranged from −0.016 to 0.057 (Table 3). Comparisons revealed that only 2 populations, Nevada and northern Utah, were significantly differentiated (FST = 0.057, P = 0.006). A 2nd population comparison was performed excluding locus VeCr07. Additional populations were significantly differentiated with the exclusion of this locus (Appendix Table 7), including Nevada and northern Utah (FST = 0.051, P = 0.008), South Dakota and northern Utah (FST = 0.034, P = 0.018), and Colorado and northern Utah (FST = 0.029, P = 0.029). South Dakota and Arizona approached significant differentiation, with FST = 0.019 (P = 0.052). The correlation test performed in CoDiDi found no significant negative correlation between HS and GST (rGH = 0.017, P = 0.50), suggesting that the population differentiation statistics reported here are accurate and not being underestimated by mutation rates. Plotting microsatellite Slatkin linearized FST values against geographic distances revealed no significant relationship (r = 0.003, P = 0.32; Appendix Figure 3). Excluding locus VeCr07 from the IBD test revealed similar results (r = 0.005, P = 0.34).
TABLE 3.
Population pairwise FST values for 7 microsatellite loci of Virginia's Warbler are shown below the diagonal. Significant pairwise comparisons are indicated in bold font and with an asterisk (P < 0.05). The estimated number of migrants (NM) between populations of Virginia's Warbler in the U.S. is shown above the diagonal.
An AMOVA on reassigned groups based on mountain range barriers revealed lower variance among individuals (20%) than within individuals (79%), with 1% explained by variance among populations (Appendix Table 8). Estimating population differentiation revealed a significant overall FST of 0.011 (P = 0.01). However, population pairwise FST values were effectively zero and none were significant, suggesting that grouping individual Virginia's Warblers according to state or geographic barrier has no tangible effect on population structure.
The best clustering determined using CLUMPAK was K = 1 according to the maximum natural log (ln) probability [lnP(−1671.2)] or K = 2 according to maximum ΔK [lnP(−1699.4); Appendix Figure 4]. One limitation of ΔK is its failure to determine when the true number of clusters is K = 1. Bar plots revealed no apparent structure; therefore, all individuals could be assigned to a single genetic cluster (Appendix Figure 5). The 2nd analysis, performed with a priori population location and correlated allele frequencies, revealed K = 1 (maximum ln probability; lnP(−1671.6)) or K = 7 (maximum ΔK; lnP(−1674.4); Appendix Figure 6). The bar plots from the 2nd analysis were indistinguishable from those of the 1st analysis, providing confidence that all individuals were clustered in 1 population.
Gene flow estimates (NM) ranged from 4.8 to 22.3, with no clear pattern of differential gene flow (Table 3). However, 2 estimates stood out. The highest effective number of migrants (NM = 22.3) was between the southern Utah and South Dakota populations, while the lowest estimated number of migrants was between northern Utah and Nevada (NM = 4.8).
Demographic Analysis
The median-joining haplotype network revealed a starburst pattern with one central haplotype shared by all populations (Haplotype 1; Figure 2). Few mutational steps separated the many less-common haplotypes surrounding the central haplotype. The starburst pattern indicates limited geographic structure and suggests a signature of rapid population expansion. Values for Tajima's D were negative for all populations (Table 1) and were statistically significant in the Nevada, New Mexico, and Utah populations (P < 0.05). The negative values indicate more nucleotide variants than expected under a model of neutrality, supporting demographic expansion. Likewise, Fu's FS test values were negative for all populations (Table 1). Globally, the hypothesis of population stasis was rejected by statistically significant negative Tajima's D (−2.47, P < 0.01) and large negative Fu's FS (−43.81) values. The estimate of tau (τ) was 0.89, indicating that the time elapsed since expansion was ~34,832 yr. The Bayesian Skyline analysis was consistent with the Tajima's D and Fu's FS results. The Skyline plot suggests that the population has been increasing and that the period of the largest population increase is concordant with the time of expansion inferred using the tau statistic (Appendix Figure 7).
FIGURE 2.
Virginia's Warbler (Oreothlypis virginiae) mtDNA haplotype network: Median-joining network of all 1,041 base pairs (bp) of the ND2 region of the mitochondrial genome. Circle sizes represent the number of individuals possessing that haplotype. Circle color and/or pattern indicate assigned population. A small, white circle without a sizing guide represents an inferred or unsampled haplotype.
Haplotype 1 was the most common haplotype, being identified in all populations and in 77 individuals (Appendix Table 9). Haplotype 4 was found in a total of 6 individuals from the Colorado, Nevada, New Mexico, Arizona, and southern Utah populations. Similarly, haplotype 7 was found in 6 individuals and was shared between Nevada and South Dakota. The remaining haplotypes were found in 1–3 individuals. The South Dakota population had 1 unique haplotype (Haplotype 30), which was found in only 1 individual.
With the exception of the Arizona population, we found no evidence for a reduction in population size signifying a bottleneck effect (P < 0.05). The 2-tailed Wilcoxon test performed for all populations did reveal a possible bottleneck in the Arizona population (P = 0.008). By excluding locus VeCr07 from the bottleneck analysis in a 2nd analysis, a 2nd population (Nevada, P = 0.03) was identified as having undergone a reduction in population size, in addition to the Arizona population (P = 0.02). However, the allele frequency distribution test for both analyses (with and without VeCr07) revealed a normal L-shaped distribution for all populations, including Arizona and Nevada, as expected under mutation-drift equilibrium.
DISCUSSION
Evidence of Continuous Gene Flow
This study investigated patterns of genetic variation in the Virginia's Warbler. To date, little has been reported on the genetics of Virginia's Warblers, making this study of increasing importance as the species continues to expand its range northward. Sakai et al. (2001) suggested that maintaining genetic diversity in a recently established population is essential for successful establishment and growth of that population. Previous studies have found that new populations in recently colonized areas on the periphery of source populations exhibit fine-scale genetic structure, lower genetic diversity, and allelic frequency gradients in the direction of range expansion (e.g., House Finch [Carpodacus mexicanus]; Hawley et al. 2006; tropical house gecko [Hemidactylus mabouia]; Short and Petren 2011; southern flying squirrel [Glaucomys volans]; Garroway et al. 2011; speckled wood butterfly [Pararge aegeria]; Harper 2011; spur-thighed tortoise [Testudo graeca]; Graciá et al. 2013; and dainty damselfly [Coenagrion scitulum]; Swaegers et al. 2014). Of the empirical studies referenced here, 4 of the 6 investigated genetic diversity and population structure in a species that underwent range expansion beginning in the 1990s, the same time period when the Virginia's Warbler was documented as having expanded its range into the Black Hills of South Dakota. Taking the generation time of Virginia's Warbler into account, our study sampled individuals from at least 7 generations after the most recent expansion. The number of generations presented here may seem low and may raise questions regarding our ability to detect the signature of a bottleneck or other effects of genetic drift. However, results from previous studies suggest otherwise (e.g., Garroway et al. 2011, Harper 2011, Short and Petren 2011, Swaegers et al. 2014). Our study provides evidence that continuous gene flow from conspecifics may prevent the loss of genetic diversity from recently established populations at the edge of a species' range. Using both mtDNA and microsatellite markers, we found a lack of population structure in Virginia's Warbler. Further, gene flow estimates suggested population connectivity throughout the species' range. We agree with Song et al. (2013), who acknowledged in their study of Light-vented Bulbuls (Pycnonotus sinensis) that high levels of gene flow in expanding populations might be attributed to strong dispersal capabilities in certain avian species.
Our data support the hypothesis that alleles from immigrating conspecifics continue to supplement the recently founded population of Virginia's Warbler in the Black Hills of South Dakota. The mtDNA haplotype and nucleotide diversities were uniform throughout all populations (Table 1). In the microsatellite data, expected and observed heterozygosities were similar among all 7 populations (Table 2). These homogenous diversity and heterozygosity estimates suggest that either the initial colonizing populations in the Black Hills and other peripheral areas were large or that consistently high levels of gene flow have maintained genetic diversity throughout the range of the Virginia's Warbler. Given the natural history and preferred habitat of the species, and the number of breeding pairs currently supported in the Black Hills of South Dakota, it seems unlikely that there was a founder event with a cessation of gene flow. However, we acknowledge that the genetic markers used may not have been powerful enough to differentiate between all possible demographic scenarios.
The range of diversity and heterozygosity estimates found in Virginia's Warbler populations are comparable with those observed in other North American birds that experienced postglacial population expansion. For instance, Milá et al. (2000) described similar mtDNA haplotype and nucleotide diversities for northern populations of MacGillivray's Warbler (Geothlypis tolmiei), concluding that postglacial expansion was responsible for the lack of variation and structure seen in that warbler. Similarly, in a recent study of Clark's Nutcrackers (Nucifraga columbiana), Dohms and Burg (2013) reported high levels of haplotype, nucleotide, and allelic diversity throughout much of the nutcracker's range, and attributed these findings and limited differentiation to population connectivity.
In expanding populations, allelic frequency gradients are often present along the axis of range expansion, as Garroway et al. (2011) found in the case of southern flying squirrels. Along with the presence of these gradients, a reduction in the effective number of alleles, or allelic richness, starting from the source population and emanating out in the direction of range expansion is sometimes observed. This phenomenon stems from the occurrence of serial founder events that are often associated with population expansion (Hewitt 1996). Our data showed no evidence of a reduction in population size throughout the axis of range expansion of Virginia's Warblers, supporting our hypothesis that a bottleneck signature would not be detected in the Black Hills population. However, the Arizona population was found to have undergone a bottleneck event. Located in the source area, this bottleneck signature could reflect a rapid historical expansion from a single refugium in a population with a low effective population size. Despite the bottleneck, allele frequencies and the number of effective alleles were comparable across all microsatellite loci and among all populations.
While significant gene flow provides the best explanation for the overall lack of structure and differentiation among populations of Virginia's Warbler, we did observe significant pairwise FST values between populations from Nevada and northern Utah. Two hypotheses may explain this result. First, the mountain ranges of the Great Basin may be acting as ‘sky islands,' with uninhabitable areas in the intermountain regions limiting gene flow between Utah and adjacent Nevada, as has been observed in other species (for example, Mountain Chickadee [Poecile gambeli]; Spellman et al. 2007). Second, samples representing the Utah geographic region were collected between the 1930s and the 1960s and may not accurately represent contemporary processes. However, a low gene flow estimate for the effective number of migrants between Nevada and northern Utah (Table 3) and the fact that we did not observe significant deviations in diversity from the other sampled populations (as one might expect between contemporary and historical samples) provide support for the sky island hypothesis. When locus VeCr07 was excluded from analyses, additional populations became significantly differentiated. While this may be indicative of genetic drift affecting the northern Utah region, we believe that this is due to loss of the information provided by the VeCr07 locus. Providing additional microsatellite markers or using a next-generation sequencing approach may help to clarify the processes that are contributing to observed patterns of connectivity in Virginia's Warbler populations.
Patterns of Rapid and Recent Range Expansion
As expected, our data suggest that the Virginia's Warbler has undergone a rapid and recent range expansion. The initial expansion dates back to ~34,832 yr ago, at which time the effects of the climatic oscillations and glacial ice sheets of the Late Wisconsin glacial stage were affecting North America (Dawson 1992). The genetic structure of species that underwent historical range expansion following the retreat of these ice sheets is characterized by a starburst phylogeny (e.g., Hairy Woodpecker [Picoides villosus]; Klicka et al. 2011; Downy Woodpecker [P. pubescens]; Pulgarín-R and Burg 2012; and Clark's Nutcracker; Dohms and Burg 2013). Like those species, Virginia's Warbler exhibit one widespread haplotype common throughout all geographic regions, with the occurrence of many singleton haplotypes emanating from that central haplotype (Haplotype 1; Figure 2). Results from Tajima's D and Fu's FS neutrality tests of mtDNA sequences further support demographic expansion in Virginia's Warbler.
One haplotype is unique to the population in the Black Hills of South Dakota. Given the recent colonization of the Black Hills and the lack of genetic structure observed across all locations, this unique haplotype could possibly be shared with the closest breeding population in the Laramie Mountains in Wyoming, which was not sampled. We do not believe that this lack of sampling alters our results, as all summary statistics were quite similar for all populations. We presume that the inclusion of individuals from Wyoming would not skew the homogeneity of values. Further, it would be difficult to determine the directionality of gene flow between the Black Hills and Wyoming populations, given the low amount of divergence in our dataset.
The data that support range expansion in the Virginia's Warbler are consistent with a natural, postglacial expansion initiated as a response to altering habitats and changing climatic conditions. In light of changing environmental conditions resulting from global climate change and other anthropogenic effects, the results from this study provide further evidence that, with species expanding or shifting their distributions, the loss of genetic diversity can be diminished by continuous gene flow. This study may also have implications for other avian species with ecological requirements similar to those of the Virginia's Warbler. For instance, Bushtits (Psaltriparus minimus) and Gray Vireos (Vireo vicinior) are 2 species that occupy similar niches and could follow the Virginia's Warbler's lead in eventually expanding their range northward. However, given that these species exhibit different dispersal capabilities we would not expect the same level of population connectivity observed here.
ACKNOWLEDGMENTS
We thank South Dakota Game, Fish and Parks and the U.S. Fish and Wildlife Service for collecting permits. We also thank the Museum of Vertebrate Zoology at Berkeley, the Museum of Southwest Biology, the Burke Museum, Louisiana State University Museum of Zoology, the Denver Museum of Nature and Science, the Utah Museum of Natural History, and the American Museum of Natural History for sending tissue samples. We are grateful to Kyle Kennedy and Mallory Ballinger for fieldwork assistance, Forrest Cain and Oxana Gorbatenko for laboratory assistance, and Dr. Brian Smith for providing edits and comments on the manuscript. A special thank you to Julian Dupuis for assistance with figure construction. The use of any trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government. The findings and conclusions in this article are those of the authors and do not necessarily represent the views of the U.S. Fish and Wildlife Service, South Dakota Fish and Game, or Black Hills State University.
Funding statement: Our field and laboratory research was partially funded by Black Hills State University, South Dakota Biomedical Research Infrastructure Network (BRIN), and the National Science Foundation (DEB0814841). The funders did not provide any input into the contents of this study. The funders did not require approval of the manuscript prior to submission and publication.
Ethics statement: All trapping, handling, and collection procedures were carried out in accordance with Black Hills State University Institutional Animal Care and Use Committee protocol A-12-06, South Dakota Fish and Game permit license #13-40, and U.S. Fish and Wildlife Service permit MB101618-0.
Author Contributions: C.M.B. and G.M.S. conceived and designed the study; C.M.B. and G.M.S. performed the experiments; G.M.S. contributed materials and resources; C.M.B. and G.M.S. analyzed the data; and C.M.B. and G.M.S. wrote the paper.
LITERATURE CITED
APPENDIX FIGURE 3.
Mantel test of isolation by distance for Slatkin linearized FST values for (A) mtDNA (r = −0.004, P = 0.39), and (B) microsatellites (r = 0.003, P = 0.32) for Virginia's Warbler in the U.S.
APPENDIX FIGURE 4.
(A) Optimal number of genetic groups (K; Evanno et al. 2005) and (B) probability of K of Virginia's Warbler in the U.S. using median values of ln P(data) as evaluated in CLUMPAK (Kopelman et al. 2015). K = 1 is suggested from the probability of K. Results are based upon 20 replicate runs per K with no a priori population information as performed in program Structure (Pritchard et al. 2000).
APPENDIX FIGURE 5.
CLUMPAK (Kopelman et al. 2015) bar plot population assignment for the number of genetic groups K = 1–5 assuming no a priori population information from 7 populations of Virginia's Warbler in the U.S.: Nevada (NV), South Dakota (SD), New Mexico (NM), Arizona (AZ), Colorado (CO), northern Utah (N. UT), and southern Utah (S. UT). Bar plots for K = 6–10 are not included. No apparent geographic structure was detected and all individuals could be assigned to 1 genetic cluster.
APPENDIX FIGURE 6.
(A) Optimal number of genetic groups (K; Evanno et al. 2005) and (B) probability of K of Virginia's Warbler in the U.S. using median values of ln P(data) as evaluated in CLUMPAK (Kopelman et al. 2015). K = 1 is suggested from the probability of K. Results are based upon 20 replicate runs per K with the assumption of correlated allele frequencies among populations and population location priors as performed in program Structure (Pritchard et al. 2000).
APPENDIX FIGURE 7.
Bayesian Skyline plot suggesting that the overall population of the Virginia's Warbler has been increasing. The thick black line indicates the effective population size through time, and the gray lines represent the 95% CI.
APPENDIX TABLE 4.
Primers used for mtDNA ND2 sequencing of Virginia's Warbler in the U.S.
APPENDIX TABLE 5.
Polymorphic microsatellite loci and primers used for genotyping Virginia's Warbler in the U.S.
APPENDIX TABLE 6.
Population pairwise genetic distances GST (upper diagonal) and FST values for mtDNA (below diagonal) of Virginia's Warbler in the U.S. There were no significant pairwise comparisons at P < 0.05.
APPENDIX TABLE 7.
Population pairwise FST values for 6 of the 7 microsatellite loci analyzed (excluding locus VeCr07) for Virginia's Warbler in the U.S. VeCr07 was omitted here for the possible presence of null alleles or allelic dropout. The remaining 5 microsatellites initially screened (Appendix Table 5) failed to amplify and were not retained for analysis. Significant pairwise comparisons (P < 0.05) are highlighted in bold font.
APPENDIX TABLE 8.
AMOVA results from reassigning individual Virginia's Warbler into 6 groups (Great Basin, Wasatch Range, eastern Rocky Mountains, western Rocky Mountains, mountains in southern Arizona and Mexico, and the Black Hills) according to mountainous areas of the U.S. Results corroborate other analyses suggesting no geographic structure in the Virginia's Warbler.
APPENDIX TABLE 9.
Geographic distribution of haplotypes throughout Virginia's Warbler sampling regions in the U.S. and Mexico.
