Grebes (Aves: Podicipedidae) are a threatened family containing species that vary widely in demography. Podicipedidae includes several species that are either rare and confined to a single lake basin, or abundantly distributed across several continents. The most speciose genus, Podiceps, particularly the eared grebe lineage, best reflects this curious demographic pattern by representing the most abundant of extant grebes, several critically endangered species, and a recently extinct species. Here, we obtained genetic data from 3 mitochondrial markers to make phylogenetic and population genetic inferences about the eared grebe clade. Using DNA from tissue, feather, skin, and toe pads, our sampling encompassed all species and subspecies, including the extinct Colombian Grebe (Podiceps andinus) and migratory and resident populations of the North American Black-necked Grebe (P. nigricollis californicus). Bayesian inference yielded novel insights into the dynamics of this group, particularly the recent ecological isolation and incipient speciation of the Colombian and Junin (P. taczanowskii) grebes, as determined from limited genetic divergence and rapid evolution of plumage color and bill shape (elongation, deepening of the culmen). DNA barcode and cytochrome b distances supported these inferences. Population genetic and divergence time analyses further revealed that the abundance of the North American Black-necked Grebe is likely associated with mid-Pleistocene dispersal from South America followed by late Pleistocene expansion during a time when hypersaline lake habitat accommodated large populations of staging birds. In conclusion, the demographic variation among species in the eared grebe group can be explained by recent ecological speciation of both a sympatric and an allopatric nature. Future investigation is warranted to determine whether this pattern of speciation and associated rapid phenotypic divergence can be extended to other grebe taxa.
INTRODUCTION
Grebes (Aves: Podicipedidae) are a small, cosmopolitan family of aquatic birds composed of a modest 22 species in 6–7 genera (Fjeldså 2004). While representative grebe taxa have been included in ordinal-level molecular phylogenetic studies, population genetic investigations have been limited to relatively small populations of the Great Crested Grebe (Podiceps cristatus) in New Zealand (Robertson and Gemmell 2002) and the Horned Grebe (P. auritus) in eastern Canada (Boulet et al. 2005). The relative lack of genetic investigation of grebes at the population level is surprising, given that grebes display extensive variation in population structure, geographic distribution, demography, and conservation status (O'Donnel and Fjeldså 1997, Fjeldså 2004). Differences in population structure and demographic history are expected to leave distinct genetic footprints. Alongside confirmation (or rejection) of expected genetic footprints, e.g., reduced genetic diversity and gene flow in demographically declining species, population genetic studies inform the tempo of genetic divergence among and within species. This temporal information is especially useful for grebes, given the paucity of reliable fossils constraining divergence time among species, and it also allows for interpretation of the rate of bill shape evolution. Bill shape in grebes is known to correspond closely with feeding specialization, and is rapidly modulated by competitive interactions (Fjeldså 1983, 2004).
Podiceps is the largest and most recently diverged genus within Podicipedidae and is often divided into 2 clades, the horned or crested grebe clade (3 species) and the eared grebe clade (5 species; Fjeldså 2004). Species in the horned grebe clade maintain relatively widespread, Holarctic distributions, whereas species in the eared grebe clade largely inhabit the New World and are marked by characteristic ear plumes on the sides of the head. Grebes in the eared clade represent a curious pattern of demography, being composed of either abundant, widespread species (including the most abundant of all grebe species, the Black-necked Grebe [P. nigricollis]), or of rare species with confined distributions. The latter category includes several critically endangered species (Fjeldså 1984), such as the Hooded Grebe (P. gallardoi) in the southern Andes that displays a bill specialized for feeding on larger invertebrates; the Junin Grebe (P. taczanowskii), which is confined to Lake Junín in Peru and is thought to have evolved from a population of Silvery Grebe (P. occipitalis) trapped in a late Pleistocene glacial refuge, leading to divergence in feeding specialization through character displacement in the bill (Fjeldså 1983); and the recently extinct Colombian Grebe (P. andinus). P. andinus was discovered in the Bogotá and Ubaté wetlands of the Eastern Andean Cordillera of Colombia in the 1920s. In 1945, this grebe was observed year-round on Lake Tota, slightly farther north in Colombia, and by the late 1950s was classified as a new subspecies of P. nigricollis. Prior to its extinction, P. andinus was classified as a full species based on differences in bill length and plumage, having chestnut coloration on its neck and ear-plumes and a grayish crown (Collar et al. 1992, Fjeldså 1993, 2004). Specialized to feed on small invertebrate prey, the rare eared grebes depend heavily on shallow lakes with very high densities of such prey, and avoid lakes with large fish populations that may control the trophic systems of wetlands. They are therefore mainly found in isolated wetlands or in regions with many ephemeral wetlands where they can move among lake basins and breed in those which are in an early successional stage. Outside the breeding season, migratory species may congregate in huge numbers in some saline lakes with extraordinary densities of brine-shrimps or brine-flies.
Genetic investigation of the highly abundant North American subspecies of the Black-necked Grebe (P. nigricollis californicus) is of special interest due to uncertainty about the extent of gene flow among populations using different migratory routes and the historical timing of population expansion. Specifically, mark–recapture studies of banded P. n. californicus have revealed 2 primary, multistaged migration pathways: (1) from breeding sites in the western Great Basin to the Salton Sea and the Gulf of California with a multimonth stopover in Mono Lake in northeastern California, and (2) from breeding sites in the eastern Great Basin to the Salton Sea and the Gulf of California with a multimonth stopover in the Great Salt Lake, Utah (Jehl and Yochem 1986, Boyd et al. 2000). Studies suggest that ~99% of all P. n. californicus inhabit these 2 staging lakes following the breeding season (Jehl and Johansson 2002, Jehl et al. 2003). A third hypothesized pathway includes direct migration to the Gulf of Mexico from breeding sites in the easternmost portion of the range (Jehl and Yochem 1986, Banks and Clapp 1987). Additional support for these pathways is limited due to the difficulty of capturing, banding, and resighting bands on a bird that performs nocturnal migrations and spends the majority of its time in and on the water (Jehl and Yochem 1986, 1987, Jehl 1990). For these reasons, genetic investigation may better elucidate the extent of genetic differentiation (i.e. long-term genetic isolation) among individuals using the 3 migratory pathways. Nonmigratory, resident populations of P. n. californicus have also been documented. However, little is known about these populations and distribution maps vary as to where they are located. Most support for a distinct population exists in the Valley of Mexico, where individuals have been documented in the highland lakes near Veracruz, Mexico, in breeding plumage during the breeding season (Loetscher 1955, Dickerman 1969). Evidence of nesting and downy young has also been recorded (Dickerman 1969, Wilson et al. 1988), but unique morphological differences are absent (Dickerman 1969). If this geographically isolated population has been relatively stable over time, it too, like the populations using the distinct migratory pathways, may harbor unique genetic variation detectable with molecular tools.
Here, we employ molecular genetic techniques and morphometric measurements of bill dimensions to: (1) determine phylogenetic and population genetic diversity in the eared Podiceps, (2) assess presence of genetic structure among different migratory and resident populations of the North American Black-necked Grebe, (3) relate diversity in bill shape to feeding specialization (see Fjeldså 1981, 1983), and (4) explore possible connections between genetic diversity and diversity in feeding ecology.
METHODS
Sampling
To investigate the phylogenetic relationships and biogeographic history of eared Podiceps, we evaluated all species and subspecies in the eared grebe clade. Tissue, skin, and toe pad samples were acquired from museum collections or in the field, including samples from 15 of the 16 known specimens of P. andinus (2 of which were used by R. Meyer de Schauensee [1959] to describe the species). Rollandia rolland, Podiceps grisegena, and P. auritus were included as outgroup taxa.
To evaluate the demographic history and genetic diversity of P. n. californicus, we acquired skin and primary-feather samples from throughout the subspecies' range, with emphasis on the staging lakes and on acquiring adequate numbers of samples from each proposed migratory pathway (Figure 1). Primary sampling locations included Mono Lake, California, and salt retention ponds near Great Salt Lake, Utah. Skin samples were acquired from museum collections and feathers were collected by rehabilitation employees of the trona industry in Green River, Wyoming. Sequences from GenBank were included in our analysis when available (Appendix Table 1).
FIGURE 1.
Distribution and sampling locations of Podiceps nigricollis californicus (Black-necked Grebe) and P. andinus (Colombian Grebe). P. n. californicus is distributed throughout North America, represented by 3 shaded regions showing its breeding range (yellow), wintering range (blue), and putative year-round locations (orange). Hypothesized migratory pathways are indicated by colored arrows, and sampling of representative taxa from these pathways are represented by colored, filled circles, where the smallest circle represents n = 1. The open black circle represents the sampling distribution of the putative resident population in the Valley of Mexico, and the open red circle represents the sampling size and former distribution of the extinct P. andinus in South America.
Genetic Markers
For phylogenetic analysis we used the full protein-coding cytochrome b (cyt b) gene (1,143 base pairs [bp]) and the protein-coding “barcode” region of the cytochrome oxidase I (COI) gene (699 bp). We used a portion (5′ end) of the noncoding control region (367 bp) to evaluate the historical demography and genetic diversity of P. n. californicus and P. andinus. Primers for cyt b and COI were designed to target small (200–250 bp) overlapping fragments (Appendix Table 2). The control region primers were designed in 2 grebes for which whole mitochondrial genomes are published in GenBank, Tachybaptus novaehollandiae (NC_010095) and Podiceps cristatus (NC_008140). Primers were designed in the conserved glutamine tRNA gene (GluF) flanking the 5′ end of the control region, and in a conserved region of the 12S rRNA gene (12S30R) flanking the 3′ end of the control region. Primers amplified only 3 of 7 tested grebe species, P. auritus, P. grisegena, and P. major, and therefore internal primers 550R and 530F were designed based on these sequences to target shorter fragments to enhance the likelihood of amplification in the eared grebe clade. Overlapping internal primers (CR3F, CR4R, and CR5F) were designed within the GluF-550R fragment and used in this study based on an alignment of the GluF-550R fragment in P. n. nigricollis and P. n. californicus.
DNA Isolation and Genetic Marker Amplification
Samples (~25 mg) were digested in 320 μl lysis buffer ATL (Qiagen, Valencia, California, USA) and 80 μl proteinase-K (10 mg ml−1) at 56°C. Following complete digestion, DNA was isolated using a standard phenol-and-chloroform–based isolation method, and DNA pellets were dissolved in an elution buffer containing Tris-EDTA. DNA isolated from the skin, feather, and toe pad samples underwent an additional purification step using the PowerClean DNA Clean-Up Kit (MoBio Laboratories, Carlsbad, California, USA).
Amplification of target DNA sequences was carried out by PCR in a Peltier thermal cycler (MJ Research PTC-225; Bio-Rad, Hercules, California, USA). PCRs were completed with an initial 3 min at 95°C, followed by 38 cycles of 30 s at 95°C, 38 s at 50°C, and 50 s at 72°C, followed by 10 min at 72°C and subsequently at 4°C. Each 50 μl reaction contained 2.5–50.0 ng of DNA, 31.0 μl of dH2O, 5.0 μl of 10× taq reaction buffer, 5.0 μl of 20 mM MgSO4, 1.0 μl of 10 mM dNTP, 0.2 μl of 5u μl−1 Taq DNA polymerase (USDNA Biotech, Fort Worth, Texas, USA), and 2.6 μl of both the forward and reverse primer at 5 mM. Amplification of target DNA was confirmed by gel electrophoresis in a 2% agarose A gel. Target DNA was purified for sequencing using 0.34 μl exonuclease, 0.68 μl shrimp alkaline phosphatase, and 0.68 μl dH2O for each 10 μl of PCR product, and run in a thermal cycler at 37°C for 30 min, followed by 80°C for 15 min and an indefinite amount of time at 4°C.
Data Analysis
Phylogenetic inference
DNA sequencing was outsourced to Macrogen (Seoul, South Korea), and sequences were manually trimmed and cleaned, primer sequences removed, and consensus sequences aligned (and verified by eye) using Sequencher 4.8 (Gene Codes Corporation, Ann Arbor, Michigan, USA) prior to analysis. Phylogenetic analysis was carried out using maximum parsimony (MP), maximum likelihood (ML), and Bayesian inference (BI) tree-building methods. MP trees for cyt b, COI, and the concatenated dataset (cyt b + COI) were constructed in PAUP* 4.0 (Swofford 2003) using a heuristic search method and a tree-bisection-reconnection (TBR) branch-swapping algorithm with 1,000 bootstrap replicates. The control region gene tree was also constructed in PAUP*; however, this tree construction was performed using the “fast” stepwise-addition method and no branch-swapping algorithm due to the large number of taxa (n = 97). This dataset was rooted with P. occipitalis and P. taczanowskii based on results (reported herein) from the cyt b and COI analyses.
ML trees for cyt b, COI, and the concatenated dataset (cyt b + COI) also were constructed in PAUP* using a heuristic search method and a TBR branch-swapping algorithm with 1,000 bootstrap replicates. Each dataset, however, was run with different parameters based on results from jModelTest 0.1.1 (Guindon and Gascuel 2003, Posada 2008) evaluated with standard Akaike's Information Criterion (AIC; Akaike 1974), which penalizes for increasing the number of parameters in the model, taking into account not only goodness of fit but also the variance of the parameter estimates (Posada and Buckley 2004; Appendix Table 3).
BI trees were constructed in BEAST 1.5.3 (Drummond and Rambaut 2007) and its companion programs BEAUTi 1.4.8, LogCombiner 1.4.8, and TreeAnnotator 1.5.3. Each dataset was run under different models of nucleotide substitution and site heterogeneity depending on results from jModelTest (Appendix Table 3). In BEAST, however, there are only 2 nucleotide substitution models available, the General Time Reversible (GTR) model (Tavaré 1986) and the Hasegawa-Kishino-Yano (HKY) model (Hasegawa et al. 1985); therefore, the best model out of these 2 options was chosen. Each dataset used empirical base frequencies, a relaxed (uncorrelated lognormal) molecular clock model (Drummond et al. 2006) without a fixed rate of nucleotide substitution, and a Yule process speciation tree model.
Bayes factor analyses were also conducted in Tracer 1.5 ( http://tree.bio.ed.ac.uk/software/tracer/) to statistically determine whether to partition each dataset by codon position to allow for different rates of nucleotide substitution at the first, second, and third codon positions, even though Bayes factors have been criticized as favoring the most parameterized (i.e. most partitioned) model (McGuire et al. 2007). Bayes factor marginal likelihoods were estimated using the method of Newton and Raftery (1994) with the modification proposed by Suchard et al. (2001). Results were analyzed using the posterior likelihood trace with 1,000 bootstrap replicates. A 2 ln Bayes factor >10 was considered highly significant in this study (Kass and Raftery 1995). Each dataset (cyt b, COI, and concatenated [cyt b + COI]) was therefore partitioned by the 3 codon positions in these analyses (2 ln Bayes factor = 295.76, 182.43, and 458.50, respectively), although partitioning did not influence topology.
For the control region gene tree, P. nigricollis and P. andinus were analyzed with the HKY + I + G model, with empirical base frequencies, a relaxed (uncorrelated lognormal) molecular clock model (Drummond et al. 2006), and a constant size coalescent tree model. Each BI analysis for every mtDNA marker was run for an MCMC chain length of 1 × 107 with a 10% burn-in, logging parameters every 1 × 103, and was repeated twice for a total chain length of 3 × 107, to ensure that independent analyses converged on the same result. The BEAST results were analyzed in Tracer 1.5 to evaluate convergence and to ensure normal distributions of parameters and effective sample size (ESS) values >200. A consensus tree was built for each dataset from the tree files of the 3 independent runs in TreeAnnotator 1.4.8, with a 10% burn-in, targeting a maximum clade credibility tree and node target heights (ages). Consensus trees for each dataset were viewed and edited in FigTree 1.2.3 ( http://tree.bio.ed.ac.uk/software/figtree/).
Divergence timing analysis
Because well-constrained fossils (see Parham et al. 2012) are unavailable for internal calibration of the eared Podiceps clade, 2 separate approaches were used to ascertain divergence times. To estimate divergence times for the entire eared grebe clade based on the protein-coding mtDNA loci, constraints were used based on node age estimates from a dataset including samples of all grebe genera (Ogawa et al. 2008). Fossil constraints on this larger dataset included the earliest known crown (modern) grebe Thiornis sociata (8.7 myr, 65 myr; Ksepka et al. 2013) and the earliest known stem flamingo Paleolodus (31 myr, 65 myr; van Tuinen et al. 2001, Ericson et al. 2006, Torres et al. 2014), which provided age range estimates for nodes in this study's analysis. Both of these fossil lineages correspond to recommended (Parham et al. 2012) justifications for fossil calibration, while the slightly older Early Miocene Miobaptus and Late Oligocene grebe material from Kazakhstan (Kurochkin 1976) remains either formally undescribed or an uncertain member of crown Podicipedidae (Ksepka et al. 2013). Priors were therefore set at 2 nodes, the root age (6.28 myr, 14.34 myr) encompassing all Podiceps–Rollandia species, and the time to the most recent common ancestor of the eared grebe clade and the Rollandia taxon set (5.82 myr, 13.42 myr).
We used a different approach to approximate sublineage ages from the faster control region, to date the time of population expansion and historical gene flow events in P. nigricollis and P. andinus. This dataset was instead analyzed with a fixed rate of nucleotide substitution of 2.05% per million years. This rate was estimated by accounting for a faster substitution rate in the control region compared with the mutation rate for cytochrome b in BEAST (following van Tuinen et al. 2008). Because this approach is likely more approximate, we discuss how divergence time inference is affected when accounting for different scenarios of time-dependent mutation rate (Ho et al. 2007, Henn et al. 2009).
Analysis of genetic diversity
Measures of molecular diversity in P. n. californicus were obtained in Arlequin 3.1 (Excoffier et al. 2005) under the Tamura and Nei model, with a gamma alpha value of 0.012 based on jModelTest results and AIC calculations. Polymorphic sites, nucleotide diversity (π), and mean number of pairwise nucleotide differences were calculated in Arlequin. Haplotype diversity was calculated in DnaSP 5.10 (Librado and Rozas 2009). To analyze the geographic distribution of the genetic diversity in P. n. californicus, the control region gene tree was evaluated and a hierarchical analysis of molecular variance (AMOVA) was performed in Arlequin. Geographic populations were defined based on hypothesized isolation and/or restricted gene flow, including the nonmigratory, resident population in the Valley of Mexico (n = 5) and individuals from the eastern portion of the range representing the putative eastern migratory pathway (n = 5). AMOVA calculations were conducted by computing a minimum spanning network among haplotypes under a pairwise difference model, with a gamma alpha value of 0.012 obtained from jModelTest and AIC calculations. A fifth population was defined to evaluate the level of genetic differentiation in the Colombian Grebe; the calculations were performed using the same criteria, although with a gamma alpha value of 0.023. Estimates of intraspecific and nearest neighbor distances based on cytochrome oxidase I were calculated using the K2P model (Kimura-2-parameter model in MEGA 4; Tamura et al. 2007), as is standard practice in DNA barcoding (e.g., Kerr et al. 2007).
Historical demographic analyses
To evaluate the historical movements and population growth of P. n. californicus, a mismatch distribution and Tajima's D (Tajima 1989) and Fu's Fs tests of selective neutrality (Fu 1997) were completed in Arlequin under default settings, with 1,000 simulations. In addition, the Bayesian skyline plot, which calculates the effective breeding population size through time, was completed in BEAST 1.5.3 (Drummond et al. 2005) and constructed in Tracer 1.5. Based on the control region phylogeny presented herein, P. andinus was included in this analysis. Control region sequence data were analyzed under a coalescent Bayesian skyline tree model and a relaxed (uncorrelated lognormal) molecular clock (Drummond et al. 2006), with a fixed rate of 2.05% substitutions per million years. Priors on the root height, P. n. californicus and P. andinus (0.48 myr, 1.74 myr), were also specified based on the estimate obtained from the control region phylogeny. The model of nucleotide substitution used for the dataset was the HKY + G model in MEGA based on jModelTest results and AIC calculations. The dataset was analyzed for an MCMC chain length of 2 × 107 with a 10% burn-in, logging parameters every 2 × 103, and was repeated for a total chain length of 4 × 107.
Bill Measurements and Analysis
Bill measurements were obtained from Podiceps nigricollis (all subspecies), P. andinus, P. occipitalis (all subspecies), P. taczanowskii, and P. gallardoi by visiting collections in Europe (Amsterdam, Berlin, Copenhagen, Leiden, Paris, Vienna), North America (Chicago, Raleigh, San Francisco), and South America (Bogotá, Buenos Aires, Lima, Medellín). Total culmen length (from the base of the mandible), width, and depth (at the base of the culmen) were measured for a total of 240 specimens archived as museum skins using 200 mm (8 inch) digital calipers with 0.01 mm accuracy (model #1478, General Tools, New York City, New York, USA). All measurements were taken by one person (M. van Tuinen) using a single set of calipers, with the exception of measurements from grebes from Colombia. These latter measurements (n = 17) were performed by P. Pulgarin using comparable digital calipers with the same accuracy, and the measurements were verified by digital re-estimation from photographs using a scale marker and the ruler tool in Photoshop 12.0 (Adobe Systems, San Jose, California, USA). Estimation error between measurers was estimated to be <5% based on triplicate repeat measurements on a subset (n = 50) of specimens. Quantitative analysis was performed in JMP 10.0 (SAS Institute, Cary, North Carolina, USA). Measurements (in mm) were first log-transformed, and ellipsoids were constructed based on 50% data point coverage. The significance of different mean values was tested using the ANOVA Tukey test.
RESULTS
Phylogenetic Inference and Divergence Times
Complete DNA sequences from 2 protein-coding mitochondrial markers were obtained for representative individuals of each subspecies to reconstruct a phylogeny across the entire eared grebe clade. The gene trees of cyt b and COI were largely in agreement with one another, as reflected by increased nodal resolution in a combined partitioned analysis (Figure 2A). In both trees, P. gallardoi fell out as the earliest diverging lineage among the eared Podiceps. Then, sister to P. gallardoi was the clade containing 2 sister species groupings, 1 with P. occipitalis and P. taczanowskii, and 1 with P. nigricollis and P. andinus. Divergence time estimates from this partitioned analysis suggested a likely Pliocene divergence of eared Podiceps (Figure 2B), with splitting into the P. occipitalis + P. taczanowskii and P. nigricollis + P. andinus clades commencing in the early Pleistocene. P. andinus formed a monophyletic group in the cyt b tree, reciprocally monophyletic to P. n. californicus, but nested within P. nigricollis as a whole. In the barcode tree, however, P. andinus was not monophyletic. P. occipitalis, P. taczanowskii, and P. nigricollis also did not form monophyletic groups in these gene trees (analyzed separately or jointly). Thus, we chose next to sample faster-evolving mitochondrial DNA sequences from a larger number of grebes. Using this extended sampling approach in the mitochondrial control region, the extinct P. andinus was placed firmly within P. n. californicus, while displaying mostly unique haplotypes (Figure 3).
FIGURE 2.
Estimates of phylogeny and divergence times in the eared grebe clade. (A) Concatenated phylogeny (cytochrome b [cyt b] and cytochrome oxidase I [COI] gene partitions) of the eared grebe clade with estimates of divergence times (time scale in millions of years, Ma), including the 95% confidence intervals estimated using program BEAST (blue node bars). (B) Summary of divergence times for the numbered nodes and node support values.
FIGURE 3.
Population genetic analysis of the North American Black-necked Grebe (Podiceps nigricollis californicus) and Colombian Grebe (P. andinus). Bayesian inference control region (367 bp) phylogeny with maximum parsimony bootstrap values, maximum likelihood bootstrap values, and Bayesian inference posterior probabilities, respectively, in brackets at nodes. Filled blue, green, and purple circles represent Black-necked Grebe individuals from the 3 migratory pathways (see Figure 1), and the open black circles represent Black-necked Grebes from the putative resident population in the Valley of Mexico. Open red circles represent Colombian Grebe individuals. Also shown are clades II and I, highlighted in the 2 insets on the right, representing 2 clades that are placed directly basal to the North American Black-necked Grebe–Colombian Grebe grouping. Clade II includes representative samples from the Old World nominate subspecies P. n. nigricollis of Europe, Africa, and Asia; Clade I includes representative Junin Grebes (P. taczanowskii) and Silvery Grebes from subspecies P. occipitalis occipitalis (Argentina) and P. o. juninensis (Peru and Colombia). Control region data suggest a lack of geographic structure in the North American subspecies of the Black-necked Grebe and paraphyly of all 4 grebe species shown.
Genetic Diversity and Historical Demography
The mismatch distribution of control region sequence data from P. n. californicus (including P. andinus) was unimodal, with no indication of multiple, genetically distinct populations. However, the coalescence of P. andinus haplotypes was more recent than that of P. n. californicus, showing a maximum separation of 2 base pairs. When analyzed separately, both mismatch distributions showed a signal of expansion but with different modes. The control region gene tree also revealed no clear genetic structure in P. n. californicus. The lack of genetic structure was supported by AMOVA calculations when comparing individuals from the nonmigratory, resident population in Mexico with all other individuals (pairwise Fst = −0.01, P = 0.53), and when comparing individuals from the eastern portion of the subspecies range excluding (pairwise Fst = −0.01, P = 0.57) and including (pairwise Fst = 0.01, P = 0.57) the resident individuals from Mexico. Despite a lack of genetic differentiation among individuals of P. n. californicus, genetic diversity analyses revealed high levels of haplotype and nucleotide diversity in P. n. californicus and P. andinus (Appendix Table 4).
Further demographic analysis of the control region sequence data supported the hypothesis that there was a population expansion in P. n. californicus. The mismatch distribution revealed that the observed data closely emulated a simulated distribution of a population that had undergone expansion, with the greatest number of individuals differing by 3 mutations. With negative values, Tajima's test of selective neutrality (D = −1.75, P = 0.02) and Fu's test of selective neutrality (Fs = −26.37, P < 0.001) also suggested a population that had undergone recent expansion. The Bayesian inference skyline plot analyses supported this interpretation (Figure 4). When taking the phylogeographic, mismatch, and BI skyline information together, it can be summarized that P. andinus coalesces deeply within the history of P. n. californicus (Figure 4, arrow 2), close to the common ancestor of all extant P. n. californicus (presumably the first to colonize North America; Figure 4, arrow 1) and before the expansion of both of these grebe (sub)species (Figure 4, arrows 3 and 4). This temporal pattern may reflect either that the recolonization of South America from the north that yielded P. andinus took place soon after the colonization of North America, or that the founding population size was large enough to preserve ancestral polymorphism. Both scenarios stress the evolutionary significance of historic dispersal. The differing signals of expansion in North America for P. n. californicus vs. in Colombia for P. andinus likely reflect unique in situ demographic responses. The disparity in timing between these events cannot be conclusively pinpointed due to uncertainty in mutation rate (Figure 4, arrows 3 and 4 vs. arrows 3' and 4'), but approximates the late Pleistocene.
FIGURE 4.
Bayesian inference skyline plot of Podiceps nigricollis californicus. Two alternative Bayesian inference skyline analyses are shown here, labeled with 4 focal demographic events for this lineage, which includes the founding of the now extinct P. andinus. Point 1 depicts the time of initial coalescence of extant P. n. californicus; Point 2 depicts the time of divergence between P. andinus and P. n. californicus; Points 3 and 4 depict the times of demographic expansion of P. n. californicus and P. andinus, respectively. From uncertainty in mutation rate calibration due to potential mutational time-dependency, events 3 and 4 are possibly younger (3' and 4') than the ages estimated from a strict molecular clock with a single mutation rate (3 and 4). Shown in light blue is the Holocene period from 10,000 yr ago to present. Effective population size is shown multiplied by generation time (× Tau).
Bill Measurements
Bill measurements indicated extensive overlap among subspecies and 3 of 5 species in all 3 dimensions (Figure 5; Podiceps gallardoi not shown). However, these measurements also highlighted the distinctly longer bill of P. taczanowskii (the Junin Grebe), with no overlap observed between this species' measurements and those of its sympatric congener, P. occipitalis (the resident Silvery Grebe; Figure 5). Secondly, P. nigricollis subspecies showed increasing bill measurements, from smallest in P. n. gurneyi to intermediate in P. n. nigricollis and largest in P. n. californicus, with some overlap among subspecies. The measurements of P. andinus extended that trend further, showing an even longer bill, with no overlap observed with any other grebe (P < 0.001) except its closest relatives in North America (P = 0.44). Alongside a longer bill, the bill of P. andinus was significantly wider and deeper (P < 0.001) than that of P. n. californicus, thus signifying a unique position in 3-D space. In identical fashion, from measurements of the few available specimens (n = 3), the bill of P. gallardoi was significantly deeper and wider (P < 0.001) than the bills of P. nigricollis and P. occipitalis, but, unlike P. andinus, trended toward a shorter bill length.
FIGURE 5.
Alternative views of the 3-dimensional morphospace of bill shape in species and subspecies of the eared grebe clade. (A) 3-D view facing the length (x axis) and depth (y axis) dimensions. (B) 3-D view facing the length (x axis) and width (y axis) dimensions. Bill length, width, and depth measurements among species and subspecies in the eared grebe clade reveal significant overlap among the Black-necked Grebe (Podiceps nigricollis) and Silvery Grebe (P. occipitalis) species and subspecies. This figure also reveals clear divergence of both the Colombian Grebe (P. andinus) and the Junin Grebe (P. taczanowskii) in bill morphology in the trend toward a longer and deeper bill. This divergence suggests that there was strong selective pressure on the founding populations to fill alternative ecological niches to already resident grebe species, supporting a relatively recent speciation event relative to the other eared Podiceps.
DISCUSSION
Phylogenetics and Bill Shape Evolution
Molecular reconstruction of the relationships among eared Podiceps supports previous hypotheses (Fjeldså 2004) in suggesting P. gallardoi as the basal lineage sister to a well-supported clade encompassing P. occipitalis, P. taczanowskii, P. andinus, and P. nigricollis. Within this clade, however, classifications are ambiguous, as species and subspecies do not show reciprocal monophyly for currently established species. On the subspecific level, this could be a result of incomplete lineage sorting due to ancestral polymorphism or limitations of the genetic markers that we chose for this analysis. Yet on the species level we would not expect to see this polyphyly except under a model of recent and rapid speciation; P. taczanowskii, for example, does not form a monophyletic group in any of the trees generated in this study, despite the fact that field studies unambiguously confirm that it is specifically distinct from the sympatric P. occipitalis (Fjeldså 1981). Interestingly, however, control region sequence data do reveal with high support that the small and isolated P. occipitalis population in the Colombian Andes represents an unrecognized divergent lineage, and therefore further investigation may reveal a novel grebe taxon. Comparisons of bill shape are thus of special interest in this regard. Colombian P. o. juninensis bill dimensions do not differ significantly from those of other conspecific populations, being most similar to measurements from the population in Junín. Yet, estimates fall at the low end of shape variation in all 3 dimensions. More study is needed to ascertain whether a genetically isolated population indeed exists, and to what extent it differs in bill shape and plumage.
In North American birds, the minimum genetic distance between 2 species at the COI barcode locus averages 4.3% (Kerr et al. 2007). Among eared Podiceps, the average genetic distance between P. occipitalis and P. taczanowskii ranges from 0.3 to 1.1%, and between P. occipitalis and P. nigricollis from 1.3 to 2.5%. These distances support recent, rapid divergence relative to other North American birds. Previous hypotheses have proposed that P. andinus may represent a relict population of P. nigricollis, suggested by its less derived and melanic plumage (Fjeldså 2004). Based on the recent divergence of P. andinus and its derived phylogenetic placement presented herein, the classification of P. andinus as a full species would necessitate the unsubstantiated (Konter 2012) elevation of P. n. nigricollis and P. n. gurneyi to species status, thus supporting a model of rapid evolution among the eared Podiceps that is at odds with taxonomic practice. Based on the totality of the data, which shows distinctness of P. andinus in mitochondrial haplotype, bill shape, ecology, and biogeography, we therefore conclude that this taxon illustrates a case of incipient speciation in grebes, with a high degree of retained ancestral polymorphism in the mitochondrial markers.
Bill length differences can evolve rapidly in grebes through competitive character displacement (Fjeldså 1981, 1983). For example, bill length in P. occipitalis is the same throughout its geographic range except at Lake Junín, where it shares habitat with the ecologically similar P. taczanowskii and R. rolland. In this study, we revealed significant overlap in bill length, width, and depth among the subspecies of P. nigricollis and P. occipitalis; however, in P. andinus and P. taczanowskii there was a trend toward a longer and deeper bill, while in P. gallardoi the trend was toward a shorter and deeper bill (see also Fjeldså 2004). These differences suggest that the evolutionary pressure to exploit an ecological niche distinct from other grebes inhabiting the same wetland is driving divergence in bill morphology and feeding ecology.
Divergence Timing
The dynamic and varied South American landscape lends itself to diverse evolutionary processes that promote divergence in avian taxa (Smith et al. 2014). Our study reveals rapid recent divergence among the eared Podiceps, and estimates for the timing of the diversification of these species fall within the Pleistocene, supporting glacial fragmentation and land bridge formation as potential drivers of grebe dispersal and speciation. An increase in diversification rates during the last million years is hypothesized to have resulted from direct fragmentation of habitat by glaciers and altitudinal migration of vegetation from climate change (Weir and Schluter 2004, Weir 2006). The distribution of P. taczanowskii on Lake Junín, a known ice-free refugium during the ice ages, lends support to the role of Pleistocene climate change in the divergence of the eared Podiceps (Fjeldså 1981, 2004).
Population Genetics
Recent studies on the population genetics of North American birds reveal that many taxa harbor cryptic genetic structure (Spellman and Klicka 2007, Klicka et al. 2011, Smith et al. 2011, Lait et al. 2012, Ralston and Kirchman 2012, van Els et al. 2012, Miller et al. 2013). Based on 2–3 hypothesized migratory pathways and evidence of nonmigratory populations, it is conceivable that Podiceps nigricollis californicus also harbors a cryptic population structure. Here, however, multiple lines of evidence reveal a large and wide-ranging population that underwent a rapid expansion with relatively unrestricted gene flow across its range. Other North American birds that have undergone recent population expansions include Semipalmated Sandpipers (Calidris pusilla), Chestnut-backed Chickadees (Poecile rufescens; Lait et al. 2012), Red-breasted Mergansers (Mergus serrator; Pearce et al. 2009), and Downy Woodpeckers (Picoides pubescens; Pulgarin and Burg 2012), with emphasis placed on Pleistocene climate change as the causal agent of expansion. Our timing estimate for the expansion of P. n. californicus supports a Pleistocene role as well.
This expansion may have been conditional on the ability of eared Podiceps to exploit the enormous uncontested food resources in hypersaline and alkaline lakes (as seen in P. gallardoi, P. nigricollis, and P. occipitalis; see Fjeldså 2004) outside the breeding season. Because of the late-Pleistocene emergence of these habitats in western North America, it has been hypothesized that the abundance of P. n. californicus is a relatively recent occurrence (Jehl 2001). The Bayesian skyline plot places the time of expansion of P. n. californicus at ~0.05–0.48 myr. The estimates of divergence time and effective population size are somewhat questionable due to uncertainty related to time dependency of the mutation rate in the mitochondrial control region (Ho et al. 2007). However, the skyline plot supports a rapid increase in population size toward the present that is further supported by the mismatch distribution and the Tajima's D and Fu's Fs tests of selective neutrality. Confirmation is needed with nuclear markers, which are less prone to time-dependent mutation rates, to reject the hypothesis that population expansion occurred during an older emergence of suitable habitat.
Given that our data suggest a panmictic population in North America, the existence of distinct migratory pathways in P. nigricollis suggests that these pathways are relatively new or that they are not strictly defined, allowing for alternative routes to suitable wintering habitat near the Gulf of Mexico. In regard to the nonmigratory population in Mexico, evidence of nesting and downy young would suggest a loss of migratory behavior, mirroring an established trend in North American migratory songbirds (Winger et al. 2014); however, unique genetic haplotypes were not present in the population, again suggesting either that the population is relatively new and unique mutations have yet to accumulate, or that this is another example of the flexible nature of the species and that some individuals forego migration to utilize favorable breeding habitat in Mexico.
Conclusions
Phylogenetic and population genetic comparisons of P. n. californicus (North American Black-necked Grebe) with P. andinus (Colombian Grebe) reveal that extraordinary population size in P. n. californicus is tied to Pleistocene climate change and likely to the appearance of hypersaline habitat in North America. The clustering of mtDNA haplotypes of P. andinus with this subspecies suggests that P. andinus represents a recolonization of South America from North America by an ancestral group of grebes at a time early on in the diversification of P. n. californicus. The morphological (plumage) and ecological (bill shape) distinctness of P. andinus, along with mostly unique DNA barcodes (albeit paraphyletic with regard to P. nigricollis), was subsequently achieved. Thus, we conclude that the now extinct P. andinus represented a newly established lineage and incipient species among Podicipedidae. Furthermore, and consistent with a tendency for rapid speciation in grebes, P. andinus may represent one of several incipient species, as is indicated by DNA barcode data on P. taczanowskii (Junin Grebe; this study) and the Aechmophorus occidentalis–A. clarkii (Western Grebe–Clark's Grebe) complex (Kerr et al. 2007). Additional multilocus sampling will be required to confirm this pattern. In summary, historic habitat change likely explains both the high present abundance of P. n. californicus and the distinctness of P. andinus and P. taczanowskii. The rapid ability of grebes to respond functionally to new habitat, due to plasticity in bill shape and flight ability, may therefore have facilitated species divergence in the Podicipedidae.
ACKNOWLEDGMENTS
We would like to thank Joseph Jehl for valuable feedback throughout the completion of this project and for establishing contact with the Wyoming trona industry, and Julie Lutz who generously provided us with feather samples. We would also like to thank Brian Arbogast, Stuart Borrett, and Steve Emslie for comments on the manuscript, and Juan Amat and the many museum collections, which included the Peabody Museum of Natural History, Yale University (YPM), Zoological Museum of the University of Copenhagen (ZMUC), The Burke Museum of Natural History and Culture at the University of Washington (UWBM), University of Michigan Museum of Zoology (UMMZ), Delaware Museum of Natural History (DMNH), Louisiana State University Museum of Natural Science (LSU), Bell Museum of Natural History, University of Minnesota (BMNH), Texas Cooperative Wildlife Collection, Texas A&M University (TCWC), North Carolina Museum of Natural Sciences (NCMN), Institute of Natural Sciences at the National University of Colombia (ICN), Academy of Natural Sciences in Philadelphia, Pennsylvania (ANSP), Field Museum of Natural History in Chicago, Illinois (FMNH), and the American Museum of Natural History, New York, New York (AMNH) for access to specimens.
Funding statement: This research was supported by start-up funds to M.V.T. at the University of North Carolina, Wilmington. The funders did not have input into the content of the manuscript, nor require approval before submission or publication.
Ethics statement: Required museum specimen holding permits were obtained for the USA by M.V.T., and for South America by J.F. and P.R.P.
LITERATURE CITED
Appendices
APPENDIX TABLE 1.
Specimen acquisition information of the taxa utilized in our genetic study, including museum collection, accession number, collection location, and specimen type.
APPENDIX TABLE 2.
Primer sequences used in this study for the amplification of target genetic markers in the eared grebe clade, including the mitochondrial DNA cytochrome b and cytochrome c oxidase I gene regions, and the 369 bp (base pairs) portion of the control region.
APPENDIX TABLE 3.
Phylogenetic tree analysis parameters of the eared grebe clade based on results from jModelTest 0.1.1 (Guindon and Gascuel 2003, Posada 2008). Cytb = the mitochondrial DNA cytochrome b; COI = the protein-coding “barcode” region of the cytochrome oxidase I gene.
APPENDIX TABLE 4.
Comparison of measurements of genetic diversity for Podiceps nigricollis californicus and P. andinus.
