Quantifiable geographic variation in DNA or morphology is often used to gauge past and present levels of population interchange and has thus helped define taxonomic boundaries, resolve evolutionary histories, and develop effective conservation strategies to preserve evolutionary diversity. We examined rangewide patterns of genetic (mtDNA and microsatellites) and morphological diversity within the Lark Sparrow (Chondestes grammacus), focusing on diversity patterns across potentially independent Western, Central, and Eastern USA regional breeding zones. Phylogenetic analysis of mtDNA sequences did not support phylogeographic patterns associated with previously characterized subspecies, though the presence of regionally specific clades suggests incipient evolutionary diversification. Regional differentiation was evidenced by significant differences in morphological traits, significant levels of genetic differentiation, reduced estimates of migration among regions, and the characterization of 2 distinct populations through Bayesian clustering. Morphometric and genetic variation distinguished Western Lark Sparrow populations historically characterized as subspecies C. g. actitus from conspecifics across a secondary cline. By contrast, regional variation between Central and Eastern populations, encompassing subspecies C. g. strigatus and C. g. grammacus, was less pronounced and was consistent with a primary cline across the American Great Plains. Our results indicate that clinal variation among populations of long-distance migratory birds may reflect incipient evolutionary divergence, secondary contact zones, and local adaptation of populations to continuously variable environments.
INTRODUCTION
Regional diversifications form the leading edge of taxonomic splits within species and, thus, are central to the evolutionary process. In most cases, diversification begins when gene flow becomes limited, predominantly tied to physical barriers separating populations (Slatkin 1987, Avise and Walker 1998, Avise 2000). The process of regional diversification is usually initiated by changes in the genetic composition of populations, with concurrent changes in morphological or behavioral variation promoting local adaptations. When population connectivity is restricted over long evolutionary time, these changes may lead to the emergence of independent evolutionary lineages (i.e. reciprocally monophyletic clades) differentiating subspecies, and, eventually, to potential speciation events (Barton and Hewitt 1989, Wakeley 2000). Evidence is mounting that, in spite of limited boundaries to gene flow, widely distributed migratory birds may experience regionally specific selective pressures across varied breeding environments, which can lead to clinal variation rather than discrete population units based on morphological, behavioral, and genetic traits (James 1991, Webster et al. 2002).
The concept of “subspecies” has long been debated, particularly with regard to its scientific validity and methods to objectively classify intraspecific variation (Zink 2004, Winker 2010). Subspecies have commonly been defined as locally adapted populations representing independent evolutionary lineages. Delineating subspecies, however, has been a subject of debate, with one camp demanding reciprocal monophyly (Avise 2000, Zink 2004, Phillimore and Owens 2006) while the other suggests a “75% rule” for the assignment of individuals to proposed subspecies (Amadon 1949, Patten and Unitt 2002). Recently, additional debate has emerged regarding whether subspecies are “geographic populations diagnosable by one or more phenotypic traits” (Remsen 2010) or simply a “tool of art” or “convenience” that delineate arbitrary parcels of a species' natural variation (Fitzpatrick 2010). Subspecies delineation, under any definition, requires that there be natural breaks in the distribution of measurable characters that are geographically concordant (Patten 2010). Such breaks may be readily apparent among allopatric conspecifics, such as forest birds of Eastern versus Western North America (e.g., Ball and Avise 1992, Boulet and Gibbs 2006). Species distributed across broad geographic regions, on the other hand, may have less cause to differ because they encounter fewer overt barriers to genetic connectivity. In these cases, intraspecific variation may exist simply as clines driven by gradual changes in environmental factors.
Evidence that populations of migratory birds are regionally differentiated has become increasingly reported using various techniques, including phylogenetics, population genetics, multivariate morphometrics, and stable isotope analysis. For example, researchers have used such techniques to map subdivisions at various points throughout the annual migratory cycle in songbirds such as Black-throated Blue Warblers (Setophaga caerulescens; Rubenstein et al. 2002, Davis et al. 2006), Yellow Warblers (S. petechia; Boulet and Gibbs 2006), and Eurasian Blackcaps (Sylvia atricapilla; Berthold et al. 1992, Perez-Tris et al. 2004). However, these species represent primarily woodland or marshland bird species that inhabit discrete habitats and, therefore, may be more prone to become geographically differentiated (Newton 2003). By contrast, populations from a large number of widespread, migratory passerine species are commonly described as distinct subspecies, too often without fully evaluating local adaptations or evolutionary independence (Ball and Avise 1992, Zink 2004, Phillimore and Owens 2006, Fitpatrick 2010, Remsen 2010). Intraspecific variation in many of these cases may better represent examples of clinal variation. Gradual changes in phenotypic and genetic traits are commonly labeled “primary clines” and are considered to be largely driven by ecological factors associated with current environmental gradients. Alternatively, abrupt and discontinuous changes often represent “secondary clines” and likely reflect historical factors associated with secondary contact of differentiated populations (Mayr 1942, Endler 1982). In the present study, we have examined breeding populations of a wide-ranging Neotropical migrant, the Lark Sparrow (Chondestes grammacus), to assess whether regional divergence among subspecies is evident within multiple quantitative markers that vary in mode of inheritance and rate of evolution.
Lark Sparrows breed across a wide distribution range (25°N–52°N, 84°W–125°W) in North America and occupy diverse but interconnected habitats, including mesic savannas in the East, xeric scrublands in the American Southwest, and Mediterranean-like savannas in California (Martin and Parrish 2000; Figure 1). Breeding across this range requires that populations of Lark Sparrows must tolerate regionally unique stressors (e.g., extreme heat, drought, spring snow) and, potentially, adapt to different environmental conditions. Populations from across this range, however, overwinter in comparatively smaller and narrower areas of California and the American Southwest through southern Mexico (Figure 1) and so could overlap during the nonbreeding season. Regional variation within the Lark Sparrow has been previously described at the morphological level. In 1827, Swainson distinguished a Western subspecies (C. g. strigatus) from Eastern conspecifics (C. g. grammacus) on the basis of the former's lighter-pigmented head and body plumage (Martin and Parrish 2000). However, Martin and Parrish (2000) expressed doubt that plumage pigmentation alone could signify regional divergence or support subspecies distinctions. Plumage variation may instead be clinal, given that pigmentation gradually changes where subspecies ranges overlap (Pyle 1997, Martin and Parrish 2000). Oberholser (1932) later used morphometric differences to describe the westernmost populations of C. g. strigatus as an “Oregon” subspecies (C. g. actitus), though this distinction has received little evaluation and less acceptance (Miller 1941, Martin and Parrish 2000). There is precedent in avian subdivisions that parallel the subspecies distributions of C. grammacus. For example, areas of sister-taxa overlap, or “suture zones,” have been previously described at various locations in the Rocky Mountains and the Great Plains (Remington 1968, Swenson and Howard 2005). Despite this, subspecies designations within C. grammacus remain arbitrary, as geographic variation in quantitative traits (e.g., morphological or genetic) within this species has not been systematically evaluated.
FIGURE 1.
Range map of Chondestes grammacus showing subspecies-defined breeding regions (solid black = Eastern; dark gray = Central; light gray = Western) and sampling locations used in this study (CA-M, CA-T, NE-K, TX-T, IL-C, OH-L). Cross-hatching in the Eastern zone indicates low densities. Winter range is outlined by a dashed black line. The approximate range of the Rocky Mountains is highlighted. Adapted from Martin and Parrish (2000); Breeding Bird Survey data, U.S. Fish and Wildlife Service; and regional Rare Bird Alerts.
The present study builds on past observations of plumage and morphological variation within Lark Sparrows, seeking to determine whether breeding populations display evidence of regional divergence consistent with current subspecies characterizations. Phylogenetic analysis of mitochondrial DNA sequences, population genetic analyses of nuclear DNA markers, and variation in morphometric traits were used to gauge patterns of geographic variation and thus assess the validity and degree of subspecies differentiation in the Lark Sparrow. During the past few years, there has been extensive debate about the use of the term “subspecies” in avian biology and the problem of diagnosing distinct population segments (James 2010). This study explores the potential relevance of continuous (primary) and abrupt (secondary) clines in delimiting subspecies and discusses the evolutionary and conservation significance of regional variation in migratory bird species. Extending this approach to other avian species would help reveal historical and evolutionary processes, however subtle, that may shape the intraspecific geographic variation of birds into discrete population units.
METHODS
Sampling
The entire range of C. grammacus was classified a priori into 3 regional breeding zones: Eastern, Central, and Western (Figure 1), delineated on the basis of previously described subspecies ranges (Oberholser 1932, Martin and Parrish 2000). Within each zone, we sampled 2 sites for a total of 6 populations studied (Figure 1). Western sites included Dye Creek Preserve, Tehama County, California (CA-T), and Hopland Research and Extension Center, Mendocino County, California (CA-M). Central sites included Brown and Duncan ranches, Tom Green County, Texas (TX-T), and Cedar Point Biological Station, Keith County, Nebraska (NE-K). Eastern sites were located at Lost Mound field station, Jo Daviess and Carroll Counties, Illinois (IL-C), and Kitty Todd Nature Preserve, Lucas County, Ohio (OH-L). Intensive surveys of the OH-L population were conducted during the 3-yr study period (2005–2007), owing to the species' endangered status in the state. Sampling at the other locations occurred in 2005–2007 at TX-T and IL-C, in 2005–2006 at CA-T, in 2006 at CA-M, and in 2007 at NE-K.
Adult and volant juvenile Lark Sparrows were captured during mornings of the breeding period (May–July) using mist nets passively deployed around favored perching sites within known breeding territories. Captured individuals were banded with a numbered leg band and, in some populations, marked with a unique color leg-band combination (100% color banded in OH-L; 87% in IL-C; 36% in TX-T). Each individual was measured across 6 morphological characters: wing chord (WC) and tail length (TL) were measured to the nearest 0.5 mm; tarsus length (TR), beak depth at nares (BD), and exposed culmen (XC) to the nearest 0.1 mm; and total weight (WT) to the nearest 0.1 g. A blood sample (10–50 μL) was drawn from the ulnar vein and stored in 1 mL of Longmire's lysis buffer (Longmire et al. 1997) for subsequent genetic analyses. Adults were sexed by examination of the cloacal region. This technique was proven successful through genetic sexing by polymerase chain reaction (PCR; 97% match, personal observation). During the study period (2005–2007), a total of 467 Lark Sparrows were captured. Most samples were from Ohio (197 samples: 36 females [F], 43 males [M], 118 juveniles [J]), followed by Texas (149 samples: 43F, 50M, 56J), California (CA-T 38 samples: 6F, 10M, 22J; CA-M 12 samples: 4F, 5M, 3J), Illinois (37 samples: 9F, 13M, 15J), and Nebraska (34 samples: 11F, 14M, 9J).
Morphometric traits of adult specimens from the ornithological collections of the University of Michigan Museum of Zoology (Ann Arbor, Michigan, USA) and Field Museum of Natural History (Chicago, Illinois, USA) were also measured and analyzed separately from contemporary samples to improve our geographic representation and ability to identify primary and secondary clines. We measured WC, TL, BD, and XC, along with beak width at nares (BW) and nares to beak tip (N2T) to the nearest 0.1 mm. The museum specimens examined included 39 Western (14F, 25M), 38 Central (10F, 28M), 84 Eastern (32F, 50M, 2 unknown), and 25 (4F, 21M) from the Rocky Mountains (New Mexico to Idaho; Figure 1).
Morphometric Analysis
Morphometric data from live-captured and museum adult Lark Sparrows were analyzed separately using JMP version 8.0.2 (SAS Institute, Cary, North Carolina, USA). All measures were natural log transformed prior to analysis (WT of live-captured birds first linearized using a cube-root function). Following Butler et al. (2007), we calculated each individual's overall size (SIZE) as the arithmetic mean of all log-transformed measures. Transformed morphometric data were analyzed using multivariate analyses of variance (MANOVA), testing population and sex as main model effects, and the interaction effect (population × sex). Age (i.e. years since collection) was also incorporated as a main effect in the morphometric analysis of museum specimens, because preserved skins may shrink over time (Winker 1993). The significance of variation explained by specimens' age and between-sex differences was assessed using an exact F-test, whereas Wilks's lambda was used to gauge variation explained by differences between populations and population × sex interactions. The minimal number of canonical correlation axes to maximally differentiate populations was assessed from the MANOVA. The canonical scores along each axis were saved for each individual and these, plus SIZE measures, were independently tested for significant differences between populations using an analysis of variance (ANOVA) with a Tukey-Kramer HSD post hoc test.
Mitochondrial DNA Analysis
DNA was extracted from blood samples by proteinase-K digestion, sodium acetate protein precipitation, and ethanol precipitation of DNA (Sambrook et al. 1989). DNA samples were resuspended in TE buffer to a concentration of ∼50 ng μL−1 for further analysis. Primer set L16743 (Tarr 1995) and H920 (5′-GTC CGG CAA CCA TTA CAC TA-3′; lab designed) was used to amplify 823 base pairs (bp) of the mitochondrial control region (D-loop). The PCR reactions were run in 25 μL containing ∼50 ng DNA, 1.5 mM MgCl2, 0.08 mM dNTPs, 0.4 μM each primer, and 1 unit of Taq DNA polymerase. The PCR profile used was as follows: denaturation at 94°C (3 min); 10 cycles of 94°C (30 s), 56°C (30 s; −0.5°C cycle−1), and 72°C (45 s); 25 cycles of 94°C (30 s), 51°C (30 s), and 72°C (45 s); a final extension at 75°C (10 min); and held at 4°C. The PCR products were ethanol-precipitated and direct-sequenced using primers L16743, H920, and 2 nested primers (H598 [5′-TTC AAC CGA CCA CTT GTA TCT G-3′; lab designed] and L437 [Tarr 1995]) with ABI BigDye Terminator chemistry (Applied Biosystems, Foster City, California, USA). Sequences were visually checked for scoring accuracy and then aligned using BioEdit (Hall 1999) to generate a consensus sequence.
Control-region DNA was sequenced from 176 Lark Sparrows across the 6 populations. These sequences included the following numbers from each population: CA-T = 33 samples (10 males, 6 females, 17 juveniles); CA-M = 12 (5M, 4F, 3J); NE-K = 27 (14M, 11F, 2J); TX-T = 37 (19M, 18F); IL-C = 28 (13M, 9F, 6J); and OH-L = 39 (22M, 17M). To minimize representation from parent–offspring pairs or siblings, juveniles incorporated into the analysis were either associated with an unmarked mother or captured independently of an adult female or another juvenile of the same haplotype (i.e. a potential sibling).
Indices of haplotype (h) and nucleotide diversity (π) were generated for each population and species-wide using Arlequin version 3.1 (Excoffier et al. 2005). A hierarchical analysis of molecular variance (AMOVA) was performed using haplotype frequencies and sequence divergence (θ statistics). Fixation indices (FCT, FSC, FST) generated through the AMOVA indicated levels of variation partitioned among breeding regions, differences between populations within regions, and differences among individuals within populations, respectively. Significance values for each index were calculated using 1,000 random permutations of haplotypes among populations. Arlequin was also used to test whether DNA mismatch distributions from populations and regions significantly violated Rogers and Harpending's (1992) model of historical sudden population expansion. This analysis generated tau, a measure of time since expansion, and a raggedness index, useful for identifying recent bottleneck events.
Probabilities that populations or regions shared polymorphisms due to contemporary gene flow or incomplete lineage sorting were calculated by coalescent analysis in MDIV (Nielsen and Wakeley 2001). The MDIV settings included the Hasegawa-Kishino-Yano (HKY) model of evolution, with Mmax = 30 or 50, and Tmax = 1 or 2. We assumed a mutation rate of 1.42 × 10−7 substitutions site−1 year−1 extrapolated from a fossil-calibrated Tetraonidae phylogeny (Johnson et al. 2007). Outputs included (1) scaled effective population sizes, θ = 2Nef × μ, where Nef is the effective female population size and μ is the mutation rate over the entire sequence (i.e. 823 bp) per year; (2) scaled migration rate, M = Nef × m, where m is the migration rate per generation; (3) divergence time, T = t1/Nef, where t1 is years since population divergence; and (4) time to most recent common ancestor, TMRCA = t2/Nef, where t2 is years since genetic coalescence. Confidence intervals were extrapolated from distributions of θ and M, but not T, because these distributions did not approach zero.
The best-fit nucleotide mutation model for our mtDNA sequence data, the Hasegawa-Kishino-Yano plus gamma model (HKY+G), was selected from 28 potential models using Modeltest (Posada and Crandall 1998). This model was incorporated into a Bayesian phylogenetic reconstruction implemented within MrBayes version 3.0 (Huelsenbeck and Ronquist 2001). This analysis consisted of 2 Markov chain Monte Carlo simulations with 4 chains each, 1 swap was attempted per run every 10 generations, chains were heated using a temperature setting of T = 0.25, and the gamma value was allowed to vary. Trees were rooted by 2 outgroup taxa from the subfamily Emberizinae for which comparable D-loop sequences were available, Camarhynchus psittaculamito and Tiaris canora (GenBank accession nos. AF110424 and AF310072, respectively). Chains were sampled every 103 generations, resulting in 104 samples taken, of which the first 25% were excluded from the run convergence analysis. Only clades with credibility values >0.50 were treated as likely phylogenetic splits.
Microsatellite Analysis
The same 176 Lark Sparrows used for mtDNA sequencing were also used in microsatellite analysis. A total of 12 microsatellite primer sets developed for other passerine species were screened for successful PCR amplification of Lark Sparrow DNA. These included Asμ09, Asμ15, and Asμ18 (Delany et al. 2000), Cuμ02 (Gibbs et al. 1999), Dpμ16 (Dawson et al. 1997), Escμ1 (Hanotte et al. 1994), Gf01b (Rasner et al. 2004), Gf05 (Petren 1998), Maμ23 (Alderson et al. 1999), Mme1 and Mme12 (Jeffery et al. 2001), and Pdoμ3 (Griffith et al. 1999). Of these, only Mme1 could not successfully amplify Lark Sparrow DNA.
Microsatellite PCRs were conducted in 20-μL volumes containing ∼50 ng DNA, 1.5 mM MgCl2, 0.08 mM dNTPs, 0.4 μM each primer, and 1 unit of Taq DNA polymerase. The PCR profiles included a denaturation step at 94°C (3 min); 10 cycles of 94°C (30 s), either 57°C (Asμ09, Asμ15, Escμ1, Gf01b, Gf05, and Pdoμ3) or 60°C (Asμ18, Cuμ02, Dpμ16, Maμ23, and Mme12) annealing temperature (30 s; −0.5°C cycle−1), and 72°C (45 s); 25 cycles of 94°C (30 s), 51°C, or 55°C annealing temperature (30 s) and 72°C (45 s); and a final extension at 75°C (10 min). Microsatellite amplifications were then held at 4°C until transferred to an ABI377 genetic analyzer (Applied Biosystems) for DNA fragment analysis.
Initial analysis indicated that Mme12 was fixed for 1 allele across all samples and was therefore omitted from further analysis. The remaining 10 polymorphic loci were first analyzed for evidence of null alleles using Cervus version 3.0.3 (Kalinowski et al. 2007). Two loci revealed significantly high predicted frequencies of null alleles (Pdoμ3 = 0.1379; Gf01b = 0.1042) and were therefore excluded from subsequent analyses. The remaining 8 microsatellite loci were used to estimate population levels of genetic diversity (allelic richness and observed [Ho] and expected [He] heterozygosities). Heterozygote deficits in relation to Hardy-Weinberg expectations were tested using GenePop version 4.0 (Rousset 2008) with a sequential Bonferroni correction for significance thresholds. GenePop was also used to calculate allelic differentiation between paired populations using exact G-tests.
The partitioning of genetic variance among breeding regions and populations was tested using a hierarchical AMOVA within Arlequin. Fixation indices of variation among breeding regions (FCT), among populations within regions (FSC), among individuals within populations (FIS), and within individuals (FIT) were calculated as weighted averages over all loci using allelic frequencies (FST-like method). An AMOVA using a stepwise model of microsatellite evolution (RST method; Slatkin 1995) was not used, because the requisite assumption of normally distributed allelic frequencies was violated by most loci (personal observation). Significance values for F-statistics generated by the AMOVA were calculated in relation to 2 × 104 permutations of alleles among breeding zones, populations, and individuals.
Geographic patterning of microsatellite genotypes was further examined using a Bayesian clustering method implemented within Structure version 2.3.1 (Hubisz et al. 2009). We allowed for admixture of populations and informed the simulated assignment by sampling location (LocPrior model). This approach increases the power of assignment for recently isolated groups without causing spurious clustering (Hubisz et al. 2009). Following a burn-in period of 106 iterations, an additional 106 iterations were used to estimate log probabilities for each possible number of population clusters (K = 1–6). For each value of K, we conducted 10 independent runs and then calculated the mean (± SE) log probability, changes in log probability between successive K values (ΔK; Evanno et al. 2005), and individuals' assignment probabilities to each cluster. To evaluate whether the LocPrior model employed might have biased the outcome of the Structure assignments, we conducted a series of LocPrior-informed simulations whereby individuals were randomly assigned to 6 populations of sizes equivalent to the actual dataset. These simulations showed the highest log probability for K = 1, as would be expected with no bias of LocPrior-informed simulations (see Figure 2). Furthermore, we performed a series of 40 independent simulations of single-individual transplants between the clusters best supported by the Structure analysis. Evaluating these simulations using the LocPrior model revealed the probability of detecting true dispersers (or their descendants) and, thus, the ability of Structure to detect demographic independence between regional populations in this species.
FIGURE 2.
Plot of mean (± SE) log probabilities of detecting K clusters in Chondestes grammacus, calculated in the program Structure for the actual dataset (solid squares) and for randomized datasets (open squares). The potential number of clusters tested ranged from K = 1 to K = 6 and included 10 iterations each.
RESULTS
Morphometric Analysis
Among live-captured specimens, a MANOVA of log-transformed measures indicated significant differences between sexes, populations, and regional breeding zones. Males differed significantly from females (F = 26.5, P < 0.001) in having longer wing chords (means: M = 87.1 mm, F = 82.6 mm; eigenvector [λ] = 2.80), lighter weights (means: M = 28.0 g, F = 28.2 g; λ = −1.19), and deeper beaks (means: M = 7.14 mm, F = 7.02 mm; λ = 0.28). Sex × population interactions showed no significant effect (Wilks's λ = 0.863, approximate F = 1.06, P = 0.38), meaning that patterns of sexual dimorphism were consistent across populations. Morphological differences among populations were significant (Wilks's λ = 0.336, approximate F = 8.88, P < 0.001), with 3 canonical axes (cC1–3) explaining 97.5% of the residual variation (Figure 3A). Variation in cC1 was explained by all measures except tarsus length (TR) and significantly differentiated all 3 breeding zones (Figure 3A). Axis cC2 was associated with residual variation in all measures except beak depth (BD) and tarsus length (TR) and differentiated TX-T from all other populations (Figure 3A). Axis cC3 explained variation in all measures except wing cord (WC) and differentiated NE-K from all other populations, including TX-T (personal observation).
FIGURE 3.
Canonical plots from MANOVAs based on multiple log-transformed morphological measures of Chondestes grammacus (for abbreviations of populations and traits, see text). (A) Circles indicating 95% confidence intervals around contemporary population means (light gray = Western, gray = Central, black = Eastern). (B) Similar data based on museum specimens from throughout the 3 ranges and, separately, from the Rocky Mountains. Brackets opposite the canonical axes show significant groupings calculated using ANOVA and Tukey test (P < 0.05).
Overall, SIZE parameters (± SE) were significantly different among populations (ANOVA: F = 12.72, P < 0.001). Post hoc tests indicated significant differentiation between larger birds in CA-T (10.69 ± 0.02) and NE-K (10.67 ± 0.02) versus smaller birds in TX-T (10.57 ± 0.01), IL-C (10.59 ± 0.02), and OH-L (10.57 ± 0.01; all comparisons between clusters P < 0.01; within clusters P > 0.90). CA-M (10.65 ± 0.03) was not significantly different from NE-K (P = 0.99), CA-T (P = 0.91), or IL-C (P = 0.38), though a small sample size could have lowered the power to distinguish CA-M individuals, which tended to be relatively larger.
The inclusion of morphometric data from museum specimens in a separate analysis expanded our representation of populations throughout the species' range, including specimens from the Rocky Mountains area located between the Western and Central regions (Figure 1). Specimens' age (years since collection; range: 15–144 yr) did not have a significant effect (F = 0.52, P = 0.79). Consistent with the analysis of contemporary samples, the MANOVA revealed significant differences between the 3 defined breeding regions (Wilks's λ = 0.614, approximate F = 5.03, P < 0.001), with 2 canonical axes (mC1-2) explaining 88.2% of the residual variation (Figure 3B). Axis mC1 was primarily explained by variation in BD and WC, whereas mC2 primarily represented residual variation in WC, TL, and XC (Figure 3B). Western, Central, and Eastern regions differed significantly along ≥1 canonical axis in all comparisons. By contrast, Rocky Mountain specimens neither formed an independent cluster nor fell intermediate to Central and Western clusters, but instead overlapped with samples from the Central region (Figure 3B).
Mitochondrial DNA Analysis
Sequences from 179 Lark Sparrows revealed variable sites at 66 of the 823 nucleotides analyzed and comprised 102 distinct haplotypes (GenBank accession nos. HQ290164–HQ290265). Populations showed high levels of haplotype diversity (h), ranging from 0.909 in CA-M to 0.980 in NE-K (Table 1), driven primarily by relatively large numbers of unique haplotypes (n = 78). Nucleotide diversity (π) was low and similar across populations, ranging from 0.006 in CA-M to 0.009 in both NE-K and TX-T (Table 1).
TABLE 1.
Genetic diversity estimates (mtDNA control region and microsatellites; Ho = observed heterozygosity; He = expected heterozygosity) for 6 breeding populations (Pop) of the Lark Sparrow located throughout 3 regional breeding zones.
The AMOVA for the mitochondrial control region revealed modest, though not significant, variation among breeding zones (4.4%, P = 0.063), significant variation between populations within regions (2.7%, P = 0.03), and significant variation among individuals within populations (92.9%, P < 0.001; Table 2). Pairwise FST estimates between populations were consistently nonsignificant between populations of the same breeding region but were significant in 7 of 12 comparisons of populations from different breeding regions (Table 3). After a sequential Bonferroni correction (α′ = 0.0036), the pairwise FST value of CA-T versus OH-L remained significantly high (P < 0.0001), whereas 2 other pairings that also happened to include CA and Eastern populations showed near-significant differences (CA-M vs. OH-L, P = 0.0039; CA-T vs. IL-C, P = 0.0039; Table 3). When samples were pooled by region (W = Western, C = Central, and E = Eastern), all pairwise FST values between regions were highly significant (W vs. C: FST = 0.013, P = 0.009; W vs. E: FST = 0.021, P < 0.001; C vs. E: FST = 0.008, P = 0.009).
TABLE 2.
Hierarchical analysis of molecular variance (AMOVA) based on mtDNA control-region pairwise differences (θ-FST) and number of different alleles within nuclear microsatellites (FST-like values).a
TABLE 3.
Levels of pairwise population differentiation in mtDNA (above diagonal; distance-based FST) and microsatellite DNA (below diagonal; θ estimates). Significant pairwise differences (P < 0.05) are boldface; asterisk indicates values that are significant after sequential Bonferroni correction (α′ = 0.0036).
As expected, MDIV-estimated migration rates were generally highest between populations within regions, though also elevated across Central/Eastern and, to a lesser degree, NE-K/Western pairings (Table 4). Estimates of time since divergence varied across paired populations but showed low divergence estimates for populations within regions (5.0–12.3 × 103 yr) and much older estimates (>23 ×103 yr BP) for most other comparisons (Table 4). Paired evolutionarily effective female population sizes ranged from 6.7 × 104 for CA-M/TX-T to 12.9 × 104 for NE-K/IL-C, with considerably overlapping 95% confidence intervals (see Table 5) and did not indicate past population bottlenecks. At the regional level, estimates of evolutionarily effective female population sizes were consistently ∼14 × 104. Bayesian phylogenetic analysis of the mtDNA sequences did not support monophyletic groupings of previously described subspecies and indicated no apparent phylogeographic patterning (Figure 4). A large percentage of sequences clustered in nonresolved polytomies with representatives from each of the 3 breeding regions. That said, 14 of the 33 nodes with clade credibility scores >0.50 were region-specific (with 5, 4, and 5 nodes for the Western, Central, and Eastern breeding zones, respectively). Although most nodes grouped only 2 sequences, there were 2 entirely Western clades of 9 and 4 sequences and an Eastern clade of 4 sequences (Figure 4).
FIGURE 4.
Rooted phylogeny of Chondestes grammacus mtDNA control-region sequences reported in this study. Breeding-zone origin is reported by the differential shading of region-specific terminal branches and clades. Values at internal nodes indicate confidence levels derived from the Bayesian analysis.
TABLE 4.
Pairwise time since divergence (above diagonal; millennia since divergence (years × 103)) and migration rates (below diagonal; migration rate (m) per year per 105 individuals) estimated using MDIV. Each pairwise distribution's mode value and 95% confidence interval are provided. Note that the values have been scaled against larger timeframes or population sizes to provide more easily comparable mode values.
TABLE 5.
Evolutionarily effective female population sizes (Nef × 104) estimated from MDIV analysis of mtDNA sequence data.
DNA mismatch distributions for populations and regions, and species-wide, showed neither departures from a sudden expansion model nor significant raggedness indices (Table 6). Tau values within populations ranged from 3.85 in OH-L to 9.75 in NE-K, and tau values within regions ranged from 5.90 in Eastern to 9.66 in Central. Assuming a mutation rate of 0.142 mutations site−1 Ma−1 for the mtDNA hypervariable region (Johnson et al. 2007), the range of estimated times since expansion for regional populations was 25–41 × 103 yr (Table 6).
TABLE 6.
Summary output of mismatch distribution analyses for each population. Tau (τ) is a measure of time since expansion; θ0 and θ1 are measures of population size at the start and end of the expansion, respectively. The sum of squared differences (SSD) between the observed data (Obs) and the best-fit model of sudden expansion were evaluated in relation to simulated expanding populations (Sim). The raggedness index (RI) is a measure of stepwise variability within the observed data and was similarly evaluated in relation to a simulated expanding population.
Microsatellite Genotyping
Genetic diversity at microsatellite loci, in terms of both allelic richness and heterozygosity, did not vary widely across Lark Sparrow populations. Allelic richness per locus ranged from a low of 6.13 at CA-M to 10.63 at TX-T, with a mean (± SD) across all populations of 9.19 ± 6.09 and a mean across all samples of 13.75 ± 8.50 (Table 1). Means for observed and expected heterozygosities across all populations were 0.603 and 0.620, respectively (Table 1). Fewer alleles and lower levels of observed and expected heterozygosity were generally found in the Western populations, both of which showed 1 allele fixed at the Maμ23 microsatellite locus (i.e. all individuals were homozygotes). Inbreeding coefficients (FIS) calculated across all loci for each population were significantly higher than zero in CA-T (0.094; P = 0.004), NE-K (0.059; P = 0.04), and TX-T (0.057; P = 0.02; Table 1). GenePop analysis on a locus-by-locus basis showed that, after a sequential Bonferroni correction, heterozygote deficit was significant only in 1 population CA-T for 1 locus (Gf05; Fis = 0.205; P < 0.001).
Allele frequencies revealed significant differences between populations from different breeding zones (all pairwise P < 0.01 following sequential Bonferroni correction except TX-T vs. OH-L). The AMOVA results from microsatellite loci indicated lower, though significant, regional differences (1.66%; P < 0.001) compared to those reported from mtDNA sequences (4.37%; P = 0.06). Decreased levels of differentiation may potentially be due to the fourfold difference in the effective size of mitochondrial versus nuclear DNA (Nei 1987). In contrast to mtDNA, variation between populations of the same region (0.25%) was nonsignificant across these markers (P = 0.15; Table 2). Pairwise population estimates of θ (FST-like) also indicated significant differences between populations of different breeding zones but no significant differences between populations of the same zone (Table 3). After a sequential Bonferroni correction, Western populations remained significantly distinct from all other populations (P < 0.007), but only TX-T versus IL-C remained significant (P = 0.004) among Central/Eastern comparisons.
Bayesian simulations conducted using Structure indicated that mean log probability scores (± SE) were higher when K was 2 (−4,766.6 ± 1.6) or 3 (−4,777.9 ± 3.1), in relation to K values of 1 (−4,829.9 ± 0.3), 4 (−4,807.4 ± 5.5), 5 (−4,818.3 ± 7.1), and 6 (−4,817.1 ± 7.2) (see Figure 2). Estimates of ΔK (Evanno et al. 2005) likewise indicated that K = 2 represented the most likely number of clusters, with samples divided into a Western cluster and a Central/Eastern cluster (Figure 5). Although the mean log probability for K = 3 was significantly lower than K = 2 (t-test, P = 0.007), partitioning the samples into 3 groups hinted at some modest differences emerging between Central and Eastern populations (Figure 5). Structure analysis of simulated dispersal events produced assignments back to “origin” clusters exceeding P = 0.50 in 20% of Central/Eastern to Western and 7.5% of Western to Central/Eastern transplants. In fact, immigrant probabilities exceeded P = 0.33, a value never reached by any individual in the true analysis, in 17 of 40 (42.5%) Central/Eastern to Western and 8 of 40 (20.0%) Western to Central/Eastern transplants.
FIGURE 5.
Assignments of Lark Sparrow populations in the program Structure. Each individual is represented by a vertical line identifying the individual's estimated membership in a particular cluster. Upper plot shows results for two clusters, which showed the highest posterior probability. Lower plot shows results for 3 clusters, each corresponding to 1 of the regional breeding zones of Chondestes grammacus. Partial probabilities of region-assignment are indicated for each individual using shades of medium gray, light gray, and charcoal, for Western, Central, and Eastern regions, respectively.
DISCUSSION
Our results show that geographically disparate populations of Lark Sparrows differ, to varying degrees, across multiple genetic and morphometric traits and suggest the presence of both primary and secondary clines across putative boundaries of previously defined subspecies. The phylogenetic analysis of mtDNA sequences did not support the characterization of Lark Sparrow subspecies as independent evolutionary lineages. Neotropical migratory passerines generally tend to have shallow and unresolved mtDNA phylogenies on account of their relatively recent diversification (Avise and Walker 1998, McKay and Zink 2010). Many of these species, including C. grammacus, may have also retained ancient variation in neutral DNA markers because they remained relatively large following isolation, thus limiting lineage sorting through genetic drift (Edwards et al. 2005, Zink et al. 2005, Milá et al. 2007, Zink 2008). For instance, Pruett et al. (2008) showed levels of genetic diversity among different Song Sparrow subspecies (Melospiza melodia) in the Pacific coast region similar to those we found throughout the Lark Sparrow. Despite little phylogenetic differentiation among subspecies of C. grammacus, there appears to be incipient genetic divergence within the species manifested as regionally specific clades, repeated patterns of pairwise population FST significance, reduced estimates of migration among regions, and the characterization of 2 distinct clusters in the Structure analysis.
Emerging genetic structure among Lark Sparrow breeding regions was consistent across analyses of nuclear microsatellite loci and mtDNA, indicating the greatest pairwise genetic differences between breeding regions, particularly Western versus all other populations. The relatively small percentage of variation explained by FST values and the limited sampling in the number of populations surveyed may instill doubt regarding the biological relevance of these estimates (see Whitlock and McCauley 1999, Waples and Gaggiotti 2006). However, coalescent analyses in Structure likewise revealed a pronounced separation of Western from Central/Eastern breeding populations (Figure 5). Furthermore, 1 microsatellite locus (Maμ23) was fixed in the California populations, reinforcing the assertion that Western populations represent an independent unit. Although inferring demographic independence solely on the basis of genetic data has been deemed speculative (see Lowe and Allendorf 2010), the Structure analysis, coalescent estimates of widely varied evolutionarily effective population sizes within and among regions, and estimates of migration rates suggest that some demographic independence exists among Western and Central/Eastern Lark Sparrow breeding populations. This is not unexpected, given that the Rocky Mountains run between these groups. The demographic independence of Eastern and Central populations is not so clear, with less support for genetic divergence among these populations (see pairwise FST, Structure results at K = 3, and migration rate estimates).
Regional demographic independence itself has biological relevance, but in the context of the present study, we sought evidence for such diversifications representing either incipient evolutionary lineages consistent with subspecies classifications or clinal variation associated with locally adapted populations. Although the mitochondrial phylogeographic patterns did not support subspecies categorizations based on reciprocal monophyly, our examination of morphological traits showed significant regional diversification consistent with results from mitochondrial and nuclear DNA population genetic analyses. Morphological patterns support the premise that the species has developed regional local adaptations, even in the face of periodic overlap or population intermixing. This is shown by the disproportionately long wings of adult birds from populations farthest from wintering grounds (i.e. OH-L, IL-C, and NE-K) and a significant east-to-west increase in beak depth (Figures 1 and 3). The former can increase efficiency of long-distance flight (Alerstam 1991), and the latter is relevant because beak depth is correlated with bite force in emberizids (Herrel et al. 2005), which can facilitate foraging on the larger and harder grass seeds of drought-prone habitats such as California savannas (Baker 1972, van der Meij and Bout 2006). Although our study did not directly assess local adaptations, the significant morphological differences we observed parallel neutral genetic differences and, thus, could play a major role in promoting regional divergence in C. grammacus.
Regional variation across Lark Sparrow breeding populations may also be promoted by differences in migratory programs. This species winters in 2 general areas, California/Baja California and central Mexico (Figure 1), which correspond to hypothetical glacial refugia (Zink et al. 2001, Swenson and Howard 2005). Like many other Neotropical migrants, Lark Sparrows may be migrating to their ancestral refugia to overwinter (Irwin and Irwin 2005, Milá et al. 2006). Indeed, our estimates of time since divergence elucidated from MDIV simulations indicated a broad range of coalescence times across regions, covering the Late Pleistocene and Early Holocene periods (5–130 × 103 yr BP), though the most recent divergences likely occurred independently within the Western and Central/Eastern regions (about 5–12 × 103 yr BP; Table 4). Genetic and morphological divergences observed between Western and Central/Eastern populations could, therefore, have arisen within independent glacial refugia and been promoted, or at least maintained, by inherent differences in the ecological and migratory demands of each region (Johnson and Cicero 2004, Milá et al. 2007). This idea is consistent with our genetic and morphometric analyses suggesting the presence of a secondary cline along the western Rocky Mountains, which has previously been described as a “suture zone” of species overlap (Swenson and Howard 2005).
Implications for Conservation
Levels of differentiation and population independence are key elements for defining groups as species, subspecies, or interconnected populations. What constitutes an avian subspecies (Fitzpatrick 2010, Remsen 2010) or whether more than a few avian subspecies are valid (Zink 2004, Phillimore and Owens 2006) is likely to be an open-ended debate. We tend to agree with the idea that subspecies should represent incipient evolutionary lineages. However, we also emphasize that intraspecific geographic variation at both the genetic and morphological levels (independent of its taxonomic status) may have important evolutionary consequences and, as a result, should have conservation significance. Recognizing this, Crandall et al. (2000) developed a model using historical and contemporary rates of genetic and ecological exchange to gauge taxonomic relationships and define conservation units. This model provides 8 different levels of population distinctiveness, from separate species to panmictic populations, each corresponding to ≥1 combination of genetic and ecological exchangeability at both contemporary and historical timescales (Crandall et al. 2000). Under 4 scenarios, paired groups are ultimately determined to be distinct populations and could be candidate subspecies (cases 3, 5a, 5b, and 6 in Crandall et al. 2000). In the case of the Lark Sparrow, we found evidence of recent genetic distinctiveness at both mitochondrial and nuclear DNA levels, particularly between Western and Central/Eastern populations, though no substantial evidence of phylogenetic splits (i.e. geographically localized clades). If these shallow genetic differences were not paralleled by ecological adaptations, they would simply signify regional variation within an interconnected population. However, in the Lark Sparrow, there are significant patterns in wing and beak morphology that likely represent adaptive clines across breeding regions.
Under Crandall et al.'s (2000) definition, the Lark Sparrow would therefore be classified as a species composed by populations with recent ecological distinction (Case 5a). This classification supports subdividing the Lark Sparrow breeding range into 2 independent management units (sensu Moritz 1994), if not separate subspecies. The Western region represents 1 management unit and, because it is diagnosable across ≥1 phenotypic trait (Remsen 2010), it lends support to Oberholser's disputed Oregon Lark Sparrow subspecies, C. g. actitus (Oberholser 1932, Miller 1941). Between Central and Eastern populations, the modest genetic and morphological differences are more indicative of a primary cline than of a typical phylogeographic split. However, this cline between C. g. strigatus and C. g. grammacus is not without conservation relevance (Fitzpatrick 2010, Remsen 2010). Even though populations of the Eastern region may be experiencing gene flow from Central populations, the associated demographic exchange may not be sufficient to stem widespread declines of Eastern populations if Midwestern prairies and savannas continue to be fragmented and degraded. These peripheral populations are also vulnerable to genetic loss and may contain locally adapted variation important for the species' evolutionary potential (Vucetich and Waite 2003).
Morphological or demographic distinction without clear phylogeographic splits has been reported in several other mid- to long-distance Neotropical migratory bird species (e.g., Ball and Avise 1992, Bulgin et al. 2003, Zink et al. 2005, Davis et al. 2006, Milá et al. 2007, Zink 2008, McKay 2009) and, therefore, seem to be the norm rather than the exception for such species. Under the definition of subspecies as evolutionarily independent lineages, variation in migratory birds may be better described by focusing on clinal patterns, rather than trying to artificially force geographic variation into discrete boxes that could potentially mislead future evolutionary studies and conservation efforts. Evaluating the breadth and depth of geographic clines through combined evidence from morphological and genetic (as well as other) characters is therefore imperative if we are to distinguish historical from ecological factors that have shaped subspecies distributions (Endler 1982). Our results demonstrate that coarse-scale delineations of discontinuous clines among populations of long-distance migratory bird species can help resolve evolutionary histories, define taxonomic boundaries, and develop effective conservation strategies to preserve evolutionary diversity.
ACKNOWLEDGMENTS
We thank the following institutions and individuals for allowing field access to Lark Sparrow populations: Oak Openings Preserve, Metroparks of the Toledo Area; Kitty Todd Nature Preserve, The Nature Conservancy; Lost Mound Unit of Upper Mississippi National Wildlife Refuge, U.S. Fish and Wildlife Service; Cedar Point Biological Station, University of Nebraska at Lincoln; D. Brown and S. Duncan, Tom Green County, Texas; Hopland Research and Extension Center, University of California–Davis; and Dye Creek Preserve, The Nature Conservancy. The University of Michigan Museum of Zoology (Ann Arbor) and the Field Museum of Natural History (Chicago) allowed access to their Lark Sparrow collections. Birds were banded under permits of R. Dawkins, Angelo State University; D. Wenny, Illinois Natural History Survey; and M. Shieldcastle, Ohio Department of Natural Resources. We thank V. Bingman, T. Herman, E. Keller, J. Noland, E. Ross, and B. Swanson for field and laboratory assistance. R. Huber provided statistical assistance. We also thank M. Hauber, K. Scribner, and two anonymous reviewers for their constructive suggestions. This project was funded through State Wildlife Grants (FFY 2005 and FFY 2006) from the Ohio Department of Natural Resources Division of Wildlife, in-kind contributions from The Nature Conservancy, and the Department of Biological Sciences at Bowling Green State University.
