Tissue Sampling and Mitochondrial DNA Sequencing

We captured thrushes on breeding territories using mist nets and prerecorded vocalizations in New York, Vermont, Maine, Nova Scotia, Newfoundland, and Labrador in June and July, 2009–2015. Birds were sampled for blood from the brachial vein, banded, and released. Blood was preserved in the field in lysis buffer (0.1M Tris-HCl, 0.1M EDTA, 0.01M NaCl, 3% SDS) and subsequently frozen at the New York State Museum (NYSM). Additional blood samples, frozen tissues, and DNA extracts preserved on FTA Whatman cards were obtained from the American Museum of Natural History, Bird Studies Canada, the Canadian Wildlife Service, NYSM, University of Alaska Museum of the North, University of Washington Burke Museum, and the Vermont Center for Ecostudies; all additional samples were from breeding birds based on capture date. Genomic DNA was extracted using a DNeasy kit (QIAGEN, Valencia, California) following the manufacturer's protocol. Residual blood and tissue samples and DNA extracts are archived at –80°C at the NYSM. We sequenced the entire mitochondrial NADH dehydrogenase subunit 2 gene (ND2, 1,041 base pairs [bp]) from our samples using primers and PCR protocols from FitzGerald et al. (2017) for a total of 186 Bicknell's Thrushes, 77 Gray-cheeked Thrushes, 55 Veeries, and 2 Hermit Thrushes (Figure 1A,  Supplemental Material Table S1). Sequences from both strands were aligned, edited, and translated using SEQUENCHER 5.4.1 (Gene Codes Corporation, Ann Arbor, Michigan).

Single Nucleotide Polymorphic Genotyping

To generate a multi-locus dataset of nuclear DNA polymorphisms, a subset of DNA extracts was sent to the Institute of Genomic Diversity at Cornell University for genotyping-by-sequencing (Elshire et al. 2011). High-quality DNA extracts (>20 ng µL^–1) were digested with the restriction enzyme PstI and a unique barcode/adapter combination was ligated to each sample prior to pooling and amplification by PCR. Fragments were pooled with amplicons from other bird species to create 3 sequencing libraries, with 95 barcoded individuals per library, and sequenced on multiple lanes of a Hi-Seq 2000 Illumina sequencer (Illumina, San Diego, California). Each library resulted in ∼200 gigabase pair (Gbp) of raw data (Barrera-Guzmán et al. 2018). The 100 bp reads were filtered using the Universal Network Enabled Analysis Kit (UNEAK) pipeline (Lu et al. 2013) implemented in TASSEL 3.0 (Bradbury et al. 2007). Adapter dimers and sequences with an ambiguous base were removed and all reads were trimmed to 64 bp. Identical reads were merged into “tags” within each barcoded individual. Using a pairwise alignment, tags with 1 bp mismatch were retained as candidate single nucleotide polymorphic (SNPs). Candidate SNP loci that were present in fewer than 75% of individuals were removed, and then individuals possessing fewer than 50% of candidate SNPs were removed. To construct genotypes based on sequencing coverage, a maximum-likelihood method (Lynch 2009) with an updated version of the script from White et al. (2013) was implemented; this method retained SNPs with an AIC ≥4 units lower than the next best reconstructed genotype. To remove paralogs (duplicated regions within the genome), SNPs with observed heterozygosities ≥0.75 were removed (Baldassarre et al. 2014, Weir et al. 2015). Finally, SNPs present in fewer than 25% of individuals and individuals with fewer than 70% of SNPs were excluded. Following filtering, we retained a dataset of 5,633 SNPs from 72 birds, including 37 Bicknell's Thrushes from locations spanning the species' entire breeding range, 23 Gray-cheeked Thrushes representing both subspecies and including individuals from Siberia, Alaska, and eastern Canada, a probable late-migrating Gray-cheeked Thrush, and a small sample of eastern and western Veeries (n = 11) (Figure 2E,  Supplemental Material Table S1).

FIGURE 2.

Analyses of dataset with all (5,633) SNPs defined 3 lineages but revealed low levels of admixture. Bicknell's Thrush = BITH, orange, n = 33. Gray-cheeked Thrush = GCTH, purple, n = 23. Veery = VEER, green, n = 11. (A) Lineage assessment using a principal coordinates analysis. The percentage of variance explained by each axis is shown in parentheses. (B) Ancestry assignment in STRUCTURE for K = 3. Each line represents a single individual. (C) Bayesian lineage assessment in SNAPP from a subset of individuals (BITH n = 4, GCTH n = 5, VEER n = 2). (D) Tess Ad-Mixer map showing the geographic locations where potential hybridization may occur for BITH and GCTH. (E) Sample localities and sizes for BITH and GCTH, and their respective breeding ranges. Samples from Alaska (n = 1) and Siberia (n = 3), the migrant from Bellechasse, Québec, and VEER are not shown.

Species Limits, Divergence Times, and Admixture Analyses

Uncorrected pairwise nucleotide divergence for ND2 was calculated using DNASP 5.10 (Rozas et al. 2003). A hierarchical analysis of molecular variance (AMOVA) was used to assess whether a significant proportion of genetic variation could be attributed to species divisions based on 20,000 randomizations of the ND2 data set. A median-joining haplotype network was constructed in the program NETWORK 4.612 (fluxus-engineering.com). Coalescent times of ND2 haplotypes were estimated for each species using BEAST 2.1.3 (Drummond et al. 2012, Drummond and Bouckaert 2015) based on the ND2 mutation rate of 0.0145 substitutions/site/million years calibrated by Lerner et al. (2011) and using a Hermit Thrush (Catharus guttatus) sequence as an outgroup. The HKY model was determined to be the best model of nucleotide substitution for each species using Partition Finder 1.1.0 (Lanfear et al. 2012) based on Bayesian information criterion. Using the Yule speciation model and both a relaxed and strict lognormal clock, we estimated Bayesian maximum clade credibility trees from an alignment that included only one sequence from each haplotype obtained from our sample of Veeries, Bicknell's Thrushes, and Gray-cheeked Thrushes. For this analysis we constrained the monophyly of Bicknell's Thrush + Gray-cheeked Thrush following Voelker et al. (2013). An initial model was run for 10,000,000 generations using all default priors and sampling every 2,000 trees. The means for gamma and kappa parameters from the initial run were used in 3 subsequent runs of 50,000,000 generations each, sampling every 10,000 trees. All files were formatted using BEAUti 2 (Drummond et al. 2012). Convergence was assessed with Tracer 1.6.0 (Rambaut et al. 2014) and the maximum clade credibility tree was determined with Tree Annotator 2.1.2 ( http://beast.bio.ed.ac.uk/treeannotator) using a 25% burn-in. Node support and mean age estimates were tallied using Fig Tree 1.4.0 ( http://tree.bio.ed.ac.uk/software/figtree/).

To examine the history and directionality of gene flow between Bicknell's Thrush and Gray-cheeked Thrush, and to estimate the size of their populations prior to and following their split we used the Isolation-with-Migration (IMa; Hey and Nielsen 2007) software to implement a coalescent model of their mitochondrial divergence using the same HKY substitution model and the ND2 mutation rate (converted to a per-locus rate by multiplying by the number of nucleotides, 1,041 bp) as in the BEAST analyses. We performed several runs of 1 million MCMC generations to refine parameter priors and then ran the program with a burn-in of 500,000 followed by 10 million generations. From the output estimates of the population mutation parameter θ we calculated population sizes using the equation N₁ = θ /4µ.

Genetic clusters based on the SNPs were examined with a standardized principal coordinates analysis (PCoA) of the covariance of the genetic distance matrix calculated in GenAlEx (Peakall and Smouse 2006). Because distance-based analyses in GenAlEx do not handle missing data well (Peakall and Smouse 2006), we removed all singleton alleles (Linck and Battey 2019) and any individual with >40% missing SNP calls, resulting in a PCoA of 65 individuals and 5,604 SNPs. Genetic structure was examined using an admixture model and correlated allele frequencies with no prior information on species or populations groupings in STRUCTURE 2.3.4 (Pritchard et al. 2000, Falush et al. 2003); 1–6 populations were examined (3 replicates each) with 100,000 cycles after a burn-in of 10,000 cycles. Ninety-five percent confidence intervals were calculated, and results were analyzed following Evanno et al. (2005). TESS Ad-Mixer (Mitchell et al. 2013) was used to spatially visualize where admixture was most likely to occur by displaying spatial interpolations of ancestry (Q matrix outputs from STRUCTURE). A species tree was created with SNAPP (Bryant et al. 2012, Bouckaert et al. 2014); BEAUTi 2 (Drummond et al. 2012) was used to format input files for SNAPP using all default priors. Because a species tree run with all individuals and all SNPs did not converge, a species tree was run with (1) a reduced dataset composed of 44 SNPs present in all 72 individuals and (2) a reduced dataset composed of all 5,633 SNPs and 11 individuals representing the geographic spread of each species. SNAPP was run for 100,000 generations, with trees sampled every 100 generations after a pre-burn-in of 10,000. Tracer 1.6.0 (Rambaut et al. 2014) was used to assess convergence. DensiTree (Bouckaert and Heled 2014) was used to export trees; the Bicknell's Thrush/Gray-cheeked Thrush split from the Veery was designated as the root.

Within-Species Genetic Structure Analyses

Individuals were placed into regional groupings, hereafter referred to as “populations.” A population was defined as at least 5 geographically proximate individuals separated from other such groups by natural biogeographic barriers, primarily large bodies of water or low-elevation areas of unsuitable habitat between mountain ranges. Geographical proximity was delimited as within 600 km because fewer than 10% of juvenile Bicknell's Thrushes disperse more than 600 km from their natal site (Studds et al. 2012). This scheme resulted in 10 geographic populations of Bicknell's Thrush and 7 of Gray-cheeked Thrush (Figure 1A,  Supplemental Material Table S1). Bicknell's Thrush samples from Cte-Nord, Quebec (n = 3), and Gray-cheeked Thrush samples from Harbour Island, Nova Scotia (n = 2), and Siberia, Russia (n = 3), were not included in population comparisons because of small sample size. Veery samples (n = 55) were used to investigate species boundaries and estimate divergence times, but ND2 population structure was not examined for this species due to small and uneven sampling.

To detect any fine-scale population structure, we ran STRUCTURE for each species individually using the same procedure as above with no prior information on population groupings for the 5,633 SNP dataset. We also examined population structure between population groupings for all species for the 5,633 SNP dataset using Hudson's pairwise fixation indices (F_ST; Bhatia et al. 2013) with 1,000 replicates. For this analysis, we excluded Gray-cheeked Thrushes from Alaska due to small sample size, and grouped Bicknell's Thrushes from Nova Scotia and New Brunswick.

Measures of mitochondrial genetic diversity, including number of haplotypes (h), haplotype diversity (H), and nucleotide diversity (π), were calculated using the programs DNASP 5.10 (Rozas et al. 2003) and ARLEQUIN 3.5 (Excoffier and Lischer 2010). AMOVAs were used to calculate how much ND2 genetic variation could be attributed to subspecies or population divisions, based on 20,000 randomizations of the dataset. Pairwise fixation indices (F_ST) between geographic populations were calculated in ARLEQUIN with 100 permutations to assess significance. The relationship between geographic distance, measured using a nearest neighbor approach in ArcMap (ESRI), and genetic relatedness (F_ST) of populations was tested using a Mantel test (1,000 permutations) in ARLEQUIN. Mismatch distributions with Harpending's raggedness index in ARLEQUIN (1,000 permutations) and DNASP tested for recent population expansions.

Pleistocene Distribution Modeling and Overlap of Refugia

To examine the history of allopatry between Bicknell's and Gray-cheeked thrush we used the niche modeling program Maxent (Phillips et al. 2006) to create modern-day species distribution models (SDMs) for those thrush species, and then hindcast those SDMs to the LGM, ∼21,000 ybp (Mix et al. 2001, Clark et al. 2009). Occurrences throughout the breeding range for each thrush species (Bicknell's Thrush: n = 274; Gray-cheeked Thrush: n = 534), modern-day geospatial layers (19 bioclimatic layers and 3 tree species, described below), and Maxent protocols are the same as in FitzGerald (2017).

Bioclimatic layers for the LGM were obtained through the WorldClim database (Hijmans et al. 2005;  www.worldclim.org) for modern-day and for 2 different Pleistocene climate scenarios: (1) NCAR-CCSM4 (National Center for Atmospheric Research-Community Climate System Model 4; Vertenstein et al. 2010) and (2) MIROC-ESM (Model for Interdisciplinary Research on Climate-Earth System Model; Watanabe et al. 2011). Both scenarios are based on coupling atmosphere, land, ocean, and ice simulations from the LGM, but the NCAR scenario is less complex and stringent compared to the MIROC scenario.

SDMs are often based solely on climatic factors, but tree species assemblages are also important, proximate drivers of bird species distributions (Lee and Rotenberry 2005, Matthews et al. 2011). SDMs that also included layers describing the distribution of balsam fir (Abies balsamea), black spruce (Picea mariana), and paper birch (Betula papyrifera) had higher AUC values (area under the receiver operating curve) and more accurately matched breeding distributions of Catharus thrush species (FitzGerald 2017). Using the same occurrence dataset for these 3 tree species as in FitzGerald et al. (2017), we projected their distributions to the LGM and then included those LGM-projected habitat suitability layers in the thrush LGM-projected distributions.

The amount of Pleistocene-projected overlap between the species was calculated using Schoener's D in ENMTools 1.3; Schoener's D views the Maxent suitability output as proportional to species abundance (Warren et al. 2008, 2010) and is calculated by standardizing each species' suitability scores across the study region to a sum of 1, summing the absolute value of the differences between suitability scores for each bird species at each pixel, then subtracting half the sum difference from 1 (Warren 2010). A value of 0 indicates no overlap, and a value of 1 means complete overlap.

Species Limits and Admixture

One bird (NYSM zt-1480) captured at Bellechasse in southern Québec on May 31 that was originally identified as a Bicknell's Thrush but possessed a Gray-cheeked Thrush ND2 haplotype and SNP genotype was presumed to be a late migrant Gray-cheeked Thrush; we discuss this bird further below but we excluded it from our analyses. One individual in the FitzGerald et al. (2017) dataset, discussed in detail in that paper and in the Discussion below, was removed from all mitochondrial analyses because of its putative hybrid ancestry, and we were unable to genotype this individual using SNPs. To assess species limits, our final datasets included 186 Bicknell's Thrushes, 77 Gray-cheeked Thrushes, 55 Veery, and 2 Hermit Thrushes (Figure 1A,  Supplemental Material Table S1) for mitochondrial analyses, and 37 Bicknell's Thrushes, 23 Gray-cheeked Thrushes, and 11 Veeries for SNP analyses (Figure 2E,  Supplemental Material Table S1).

The average number of pairwise substitutions between Bicknell's Thrush and Gray-cheeked Thrush was 2.31% (range: 1.83–2.69%), similar to the averages between Veery and Bicknell's Thrush (1.82%, range: 1.54–2.31%) and between Veery and Gray-cheeked Thrush (average: 2.07%, range: 1.44–2.88%). The hierarchical AMOVA based on pairwise genetic distances between Gray-cheeked and Bicknell's thrush populations showed that 92% of the variance was explained by species designations, whereas only 6% of the variance was explained by within-population genetic differences (P < 0.0001; df = 2 and 17). The gene genealogy for ND2 revealed 3 major lineages corresponding to Bicknell's Thrush, Gray-cheeked Thrush, and Veery (Figure 1C). Node support values and coalescence estimates were highly congruent among the 3 randomly seeded BEAST runs ( Supplemental Material Table S2). The posterior probability values at the node splitting Veeries from the Bicknell's/Gray-cheeked Thrush clade ranged from 0.850 to 0.857. The median age of the coalescent node for 34 Gray-cheeked Thrush haplotypes was 242–243 thousand years (ky) (based on the 3 runs; 95% CI: 37–720 ky). The coalescent age of Bicknell's Thrush (n = 29) was slightly younger, with a median age of 188–192 ky (95% CI: 23–571 ky), whereas Veeries (n = 28) coalesced 243–249 ky (95% CI: 31–753 ky). The Bicknell's/Gray-cheeked Thrush clade coalesced 417–423 ky (95% CI: 59–1240 ky). Running the model with a strict lognormal clock reduced the ranges of all confidence intervals but did not substantially change the estimated mean divergence times. Estimates from IMa of gene flow between Bicknell's Thrush and Gray-cheeked Thrush following their split from a common ancestor were symmetrical and essentially 0 (high point m1 = 0.0003, 95% HPD = 0.0013–0.1482; high point m2 = 0.0003, 95% HPD = 0.0008–0.1093), translating to 10^–9 migrants per generation. Estimated values for the ancestral population mutation parameter q_a did not converge, but values of q₁ and q₂ produced estimated population size of 506,000 for Bicknell's Thrush and 957,000 for Gray-cheeked Thrush. We suspect the latter is an underestimate because we did not sample Gray-cheeked Thrush densely throughout its entire geographic range.

Multi-locus genotype analyses consistently revealed 3 genetic groupings. The PCoA segregated 3 groups (Figure 2A), with Veeries genetically distant from both Bicknell's Thrushes and Gray-cheeked Thrushes. STRUCTURE identified K = 3 as the most likely number of distinctive groups (Figure 2B,  Supplemental Material Figure S1), and the Bayesian SNAPP tree revealed 3 separate lineages, with Bicknell's Thrush and Gray-cheeked Thrush as sister taxa (Figure 2C). At K = 3, STRUCTURE revealed 3 Bicknell's Thrushes (from the populations of Bellechasse, Saguenay, and Nova Scotia) that had >2% (range: 2.8–11.0%) genetic admixture from Gray-cheeked Thrush; 2 of these had admixture 95% confidence intervals that did not overlap with 0 ( Supplemental Material Table S3). No Gray-cheeked Thrush or Veery showed >2% mixed ancestry. TESS Ad-Mixer showed that admixture was most likely to occur near the species' range boundaries, in northern Nova Scotia, southwestern Newfoundland, and in southeastern Québec, (Figure 2D), although admixed individuals were only found in the Bicknell's Thrush range.

All between-species population comparisons using Hudson's F_ST showed high differentiation (min F_ST = 0.234; Figure 3,  Supplemental Material Table S4). F_ST pairwise comparisons were higher between Veery/Bicknell's Thrush populations (range F_ST = 0.377–0.454) and Veery/Gray-cheeked Thrush populations (range F_ST = 0.328–0.387) than Bicknell's/Gray-cheeked Thrush population (range F_ST = 0.234–354). In general, Gray-cheeked Thrush from Siberia had slightly higher F_ST values (range F_ST = 0.279–0.354) than those from eastern Canada (range F_ST = 0.234–0.321) when compared to Bicknell's Thrush populations.

The presumed late migrant (NYSM zt-1480) possessed a Gray-cheeked Thrush SNP genotype and was not admixed. This bird was heavier (34.5 g), had a longer wing chord (105.0 mm), and a longer exposed culmen (14.3 mm) than most Bicknell's Thrushes (Frey et al. 2008; measurements from Y. Aubry, personal communication). Because it possessed a Gray-cheeked Thrush haplotype, non-admixed genotype, and Gray-cheeked Thrush morphological characteristics, its status as a misidentified late migrant was confirmed.

Within-Species Population Structure

A total of 63 ND2 haplotypes were detected in 263 Bicknell's Thrush and Gray-cheeked Thrush ND2 sequences (Figure 1B, Table 1), and a further 28 haplotypes were found in 55 Veery samples. Pairwise distances within species averaged 0.30% (maximum: 0.96%) for Gray-cheeked-Thrush, 0.12% (maximum: 0.48%) for Bicknell's Thrush, and 0.29% (maximum: 0.96%) for Veery. We found high levels of ND2 sequence polymorphism within our sample of Gray-cheeked Thrush and slightly lower levels within Bicknell's Thrush (Table 1). Gray-cheeked Thrushes had more haplotypes, higher haplotype diversity, a greater number of pairwise differences, and greater nucleotide diversity compared to Bicknell's Thrushes. Tajima's D and Fu's F were significantly negative for both species, indicating recent population expansions. The mismatch distribution for Bicknell's Thrush was smooth and unimodal, and the sum of squared distributions and Harpending's raggedness index were significant, further supporting a recent population expansion. The Gray-cheeked Thrush mismatch distribution was bimodal, with peaks at 2–3 and 6 pairwise differences. The mismatch sum of squared differences and Harpending's raggedness index were low and not significant. Together, these results indicate that Gray-cheeked Thrush consists of 2 populations that underwent recent population expansions.

The most haplotype diverse Bicknell's Thrush population was from the Gaspé Peninsula in Québec (Table 1), while the Nova Scotia population was the most divergent, with 8 significant F_ST values out of 9 comparisons (Figure 3,  Supplemental Material Table S5). All Bicknell's Thrush populations south of the St. Lawrence River were not significantly different from each other, except the Bellechasse, Québec, population, which showed slight difference with the S. New York (F_ST = 0.087, P < 0.05) and N. New York (F_ST = 0.052, P < 0.05) populations. A hierarchical AMOVA based on pairwise genetic distances comparing populations south (S. NY, N. NY, VT, NH, ME, and Bellechasse) and north (Saguenay, Gaspé, NB, and NS) of the St. Lawrence River (see Figure 1A) revealed that 10.0% of genetic variance was explained by these groupings (P < 0.01, df = 2 and 10). New Brunswick was not significantly different from any population, but this may be due to a small sample size (n = 11); removing New Brunswick slightly increased the amount of genetic variance explained by the north/south groupings (12.8%, P < 0.01, df = 2 and 9). An AMOVA comparing the Nova Scotia population to all other populations grouped revealed that 15% of genetic variance was explained among these groupings (P < 0.01, df = 2 and 10) compared to within populations (83%, P > 0.1).

FIGURE 3.

Pairwise population comparisons of fixation indices (F_ST) for mitochondrial ND2 (above the black diagonal) and 5,633 anonymous SNPs (below the black diagonal). F_ST values are indicated within each pixel. Gray pixels indicate no comparison. Between-species comparisons were not analyzed for ND2. Note that population inclusion and sizes were not identical for ND2 and SNP comparisons (see text). Population names can be found in Figure 1 and Table 1.

Analyses using the SNP datasets, however, did not reveal a similar north/south grouping within Bicknell's Thrush. The highest levels of genetic differentiation almost always consisted of a pairing with either the Nova Scotia/New Brunswick or the northern New York population (Figure 3,  Supplemental Material Table S5). The single-species STRUCTURE analyses of SNP data did not reveal significant groupings; however, individuals from Nova Scotia tended to segregate out when K = 3, and 2 individuals from the Gaspé Peninsula were distinctive when K = 4 ( Supplemental Material Figure S2A). We found no major genetic breaks corresponding to geography in the ND2 haplotype network (Figure 1B) or SNP ancestry plot ( Supplemental Material Figure S2A) and a Mantel test revealed large (r = 0.677) and significant (P < 0.01) isolation by distance among Bicknell's Thrush populations.

Ancestry plots revealed that Gray-cheeked Thrushes from Alaska/Siberia (C.m.aliciae) were distinct atK = 2, and individuals from western Labrador (C. m. aliciae) tended to segregate from Newfoundland/southern Labrador (C. m. minimus) when K = 3 ( Supplemental Material Figure S2B). Likelihood scores, however, did not specify a significant number of groupings. The Siberian population was moderately differentiated from the Newfoundland/southern Labrador (C. m. minimus) populations based on Hudson's F_ST (range: F_ST = 0.091–0.107; Figure 3,  Supplemental Material Table S4). Newfoundland/southern Labrador populations were similar to each other (range: F_ST = 0.029–0.036), and western Labrador showed higher similarity to Newfoundland/southern Labrador (range: F_ST = 0.025–0.031) than to Siberia (F_ST = 0.073). A Mantel test of isolation by distance for ND2 was low and not significant (r = 0.315, P > 0.1) when including all Gray-cheeked Thrush populations. Note, however, our analyses did not include any intervening populations between western Labrador and Siberia. Taken together, the nuclear and mitochondrial results are consistent with a scenario in which C. m. minimus and C. m. aliciae diverged followed by subsequent nuclear gene flow between C. m. minimus and geographically adjacent populations of C. m. aliciae.

TABLE 1.

Haplotype diversity measures for populations of Bicknell's Thrush, Gray-cheeked Thrush, and Veery based on complete sequences of ND2. Sample sizes are shown in parentheses. Tests involving statistical significance were not completed for populations with sample sizes <5. Number of haplotypes (h), number of private haplotypes (hpriv), haplotype diversity (H), average number of pairwise differences (k), nucleotide diversity (π), Tajima's D, Fu's F, the expansion parameter τ, mismatch sum of squared deviation (SSD), Harpending's raggedness index (r). Significance level is indicated by asterisks (*<0.05, **<0.01). Nfld refers to the island of Newfoundland, a subregion of the province of Newfoundland and Labrador, Canada.

Veeries had high levels of genetic diversity based on ND2 despite uneven and sparse sampling (Table 1). This is perhaps unsurprising considering the large geographic range between sampling localities and the fact that Veeries have several recognized subspecies (Heckscher et al. 2017) and 3 of these subspecies were sequenced in this study: the eastern Veery (C. fuscescens fuscescens), the western Veery (C. f. salicicola), and the Newfoundland Veery (C. f. fuliginosus). At K = 2, the inferred ancestry plots split eastern and western Veeries ( Supplemental Material Figure S3C), but Hudson's F_ST between eastern and western Veeries was fairly low (F_ST = 0.045; Figure 3,  Supplemental Material Table S5). We do not discuss Veery population structure further due to sparse sampling.

Pleistocene Refugia

Modern-day distribution models for balsam fir, black spruce, and paper birch had AUC values >0.85, indicating that each model fit the occurrence data well, so these SDMs were projected to the LGM of the Pleistocene ( Supplemental Material Table S6). These species are broadly overlapping, and each scenario hindcasted similar LGM refugia locations for all 3 tree species ( Supplemental Material Figure S3). All modern-day SDMs for thrushes had high (>0.89) AUC values, indicating a high fit with the occurrence data, and were subsequently projected to the LGM ( Supplemental Material Table S6). Both LGM scenarios show that Gray-cheeked Thrush had suitable habitat in Beringia, the central United States, and on the Grand Banks off the east coast of present-day Newfoundland (Figure 4). The CCSM4 scenario additionally identified suitable habitat near present-day New Jersey, Haida Gwaii, and north of the Northwest Territories. The Bicknell's Thrush MIROC scenario showed a large swath of suitable habitat in the southern United States. The CCSM4 scenario revealed suitable habitat on the continental shelf off the coast of Newfoundland and near present-day Maryland, Virginia, and North Carolina. Estimates of overlap at the LGM between thrush species revealed only 8.2% overlap based on the MIROC scenario and 12.2% overlap in the CSSM4 scenario. Under the MIROC scenario, the overlap occurred in small patches from present-day Oklahoma east to Tennessee, whereas overlap for the CCSM4 scenario was found on the Grand Banks off Newfoundland and western North Carolina. As reported in FitzGerald (2017), the distribution of balsam fir was by far the highest contributor to the Bicknell's Thrush models, followed by precipitation of the warmest quarter and the distribution of black spruce. Gray-cheeked Thrush LGM models were driven more or less equally by the distribution of black spruce and summer temperature variables.

Species Limits and Admixture

We investigated species boundaries, admixture, gene flow, and geographic genetic structure within the Gray-cheeked Thrush/Bicknell's Thrush species complex utilizing the largest to-date sampling of individuals and genetic loci from these species and their closest relative, the Veery. The relationships and species boundaries in this clade have been the subject of extensive study by ornithologists for over a century, most prominently in the examinations of voice, morphometrics, and plumage color by Wallace (1939), Ouellet (1993), and Marshall (2001), and in molecular phylogenetic studies of the genus Catharus by Outlaw et al. (2003), Klicka et al. (2005), Winker and Pruett (2006), and Voelker et al. (2013). With the exception of a single individual that is a putative hybrid (FitzGerald et al. 2017; discussed below), our ND2 sequences revealed that breeding individuals of the species were reciprocally monophyletic and are consistent with the Voelker et al. (2013) study which found that Gray-cheeked and Bicknell's thrushes are sister species. Our estimate of the coalescent times of the mitochondrial lineages (median = 420 kbp, mean = 518 kbp;  Supplemental Material Table S2) agrees broadly with phylogenetic divergence dates calculated by Weir and Schluter (505 kbp; 2004) and Voelker et al. (410 kbp; 2013) and places the Bicknell's/Gray-cheeked Thrush divergence in the late Pleistocene. The dominant pattern of reciprocal monophyly, very few admixed individuals, lack of historic gene flow following their divergence in the Pleistocene, and our finding that 5,633 anonymous SNP loci consistently defined 3 groupings that correspond to the 3 species, support current classifications that recognize these as species-level taxa.

STRUCTURE plots of 5,633 SNPs revealed no interspecific genetic admixture for Gray-cheeked Thrushes or Veeries, but did show low levels of admixture (2.8–11.0%)in 3 of 37 Bicknell's Thrushes. Bicknell's Thrush populations were genetically more similar to Gray-cheeked Thrushes from nearby populations than to Gray-cheeked Thrushes from Siberia (Figure 3,  Supplemental Material Table S4), and genetically admixed individuals were more likely to occur near the species' range boundaries (Figure 2D), suggesting that genetic admixture is better explained by a history of occasional hybridization events rather than shared ancestral polymorphism. The potential for occasional contact between the 2 species near their range boundaries is high because second-year Bicknell's Thrushes can disperse up to 700 km from their natal areas (Studds et al. 2012), a distance that likely far exceeds the gap between the species. In the early 20th century, the closest populations of Bicknell's Thrushes and Gray-cheeked Thrushes may have been only 60 km apart along the north shore of the Gulf of St. Lawrence (Marshall 2001), and the gap is apparently shrinking in Nova Scotia, where Gray-cheeked Thrushes appear to have recently (since the 1990s) colonized some islands along Nova Scotia's eastern coast (Whittam 2015). The discovery of Gray-cheeked Thrushes likely breeding to the south of Bicknell's Thrushes that breed in northern Nova Scotia may be a harbinger of increased secondary contact between the species.

FIGURE 4.

Models of Bicknell's Thrush (red/orange) and Gray-cheeked Thrush (purple/blue) distributions under 2 LGM scenarios, CCSM4 (red or purple) and MIROC (orange or blue). Darker colors show higher suitability. The insets show habitat suitability during the LGM for the central United States and Newfoundland. The gray shows the extent of the land during the LGM; modern political boundaries, coastlines, and regions discussed in text are shown for reference.

Perhaps more important than the potential contact between species at their range boundaries is the potential for introgression as migrating Gray-cheeked Thrushes pass through the breeding range of Bicknell's Thrush each spring. We identified a probable late-migrating Gray-cheeked Thrush captured in southern Québec that had <1% (95% CI: 0.0–0.7%) genetic admixture from Bicknell's Thrush using genotyped SNPs ( Supplemental Material Table S3), supporting the fact that some Gray-cheeked Thrushes are still migrating during the Bicknell's Thrush breeding season (Winker 2010, Townsend et al. 2015, Whitaker et al. 2018). Similarly, Marshall (2001) identified 2 other “misplaced” individuals: a male Gray-cheeked Thrush specimen (Canadian Museum of Nature 8298) collected in Percé, Gaspé, Quebec, on June 14, 1915, and a female Gray-cheeked Thrush specimen (Museum of Comparative Zoology 324734) collected on Mt. Washington, New Hampshire, on July 11, 1884. The breeding seasons of these 2 species overlap substantially; Gray-cheeked Thrushes arrive at breeding sites from mid-May to early June and have reported egg dates from June 6 to July 24 (Whitaker et al. 2018), and Bicknell's Thrushes arrive in mid- to late May and have reported egg dates from June 3 to July 28 (Townsend et al. 2015). Quay (1986) showed that thrushes can accumulate sperm in cloacal protuberances during migration before arrival on the breeding grounds, giving rise to the potential for cross-species fertilizations. The mating system of Gray-cheeked Thrush has not been studied, but the well-studied Bicknell's Thrush population that breeds in Vermont appears to be polygynandrous; both males and females mate with multiple partners and 69% of broods are sired by multiple males (Goetz et al. 2003). This propensity for extra-pair copulations may further increase the potential for occasional interbreeding as Gray-cheeked Thrushes migrate through the Bicknell's Thrush breeding range when both species are physiologically reproductive. The presence of 3 heavily backcrossed, admixed Bicknell's-like individuals within the range of Bicknell's Thrush is consistent with this scenario. In contrast, backcrossed individuals were not detected within the range of Gray-cheeked Thrush. Thus, despite the potential for geographic and temporal overlap between the species during the breeding season, hybridization apparently occurs very rarely, as evidenced by the mitochondrial and SNP analyses presented here. Increased sampling of individuals near their range limits could offer insight into the frequency of hybridization.

FitzGerald et al. (2017) identified one early generation Bicknell's/Gray-cheeked Thrush hybrid (NYSM zt-1278) based on mitochondrial and nuclear genetic evidence. This bird, captured on July 4, 2015, in southern Labrador, was phenotypically indistinguishable from Gray-cheeked Thrushes but carried an ND2 haplotype (not included in analyses in this study) that was widespread and common in Bicknell's Thrush (FitzGerald et al. 2017). Unfortunately, we failed to obtain SNP genotype data from this individual and are not able to resolve the question of its ancestry with this dataset. We performed a second set of IMa analyses of ND2 sequences in which we included the putative hybrid from southern Labrador. The estimated rate of gene flow jumped 2 orders of magnitude to m₂ = 0.032 (95% HPD = 0.0077–0.1618), translating to 10^–7 introgressions per generation, while the values for all other parameters remained essentially unchanged. Thus, even if we include the putative hybrid in our analyses we still find essentially zero gene flow between species following their split from a common ancestor, though we caution that the gene flow and population size estimates are based on the coalescence of a single locus. This may explain why our estimate of the global population size of Bicknell's Thrush (∼500,000) is substantially higher than the 120,000 estimated by Hill and Lloyd (2017) from Mountain Birdwatch survey data and a Bayesian model of abundance and density, predicted by elevation latitude and a forest cover layer.

The chance discovery of this cryptic hybrid documents occasional interbreeding and highlights the importance of continued sampling and genetic analysis of individuals in this species complex. Martinsen et al. (2018) reported the only other genetically supported case of hybridization in this species complex when they documented a Veery/Bicknell's Thrush hybrid from Vermont that resembled a Veery in plumage and carried a Veery mitochondrial DNA haplotype, but sang songs of both species and was heterozygous for intron alleles that segregated between the 2 species. The close resemblance of Catharus species and the lack of sexual dimorphism in the group might mask the frequency of hybridization relative to groups such as the wood warblers (Parulidae), for which many interspecific and even intergeneric hybrids are known (McCarthy 2006). If real, the rarity of interbreeding between Gray-cheeked and Bicknell's thrushes is noteworthy and begs the question as to what isolating mechanisms are keeping these species from hybridizing.

In addition to allopatric breeding ranges, other prezygotic isolating mechanisms may inhibit interbreeding. For example, the Bicknell's Thrush breeds in balsam fir forests whereas the Gray-cheeked Thrush tends to breed in black spruce scrub, although the Gray-cheeked Thrush can use other habitat types, including balsam fir forests (Whitaker et al. 2015, FitzGerald 2017, Ralston et al. 2019). However, no work to date has examined the potential for behavioral isolation between the Bicknell's Thrush and Gray-cheeked Thrush. These species sing songs that vary slightly but have nearly identical calls (Ouellet 1993, Marshall 2001). They also overlap morphologically and can have similar plumage coloration, depending on the time of year (Ouellet 1993, Marshall 2001, Frey et al. 2008). A monumental study of thrush behavioral isolating mechanisms found that birds recognized conspecifics based on song type rather than plumage (Dilger 1956), but because the Bicknell's Thrush and Gray-cheeked Thrush were considered conspecific at the time, the Gray-cheeked Thrush was not included in the study. Observing natural behavioral interactions between the species is rare due to the difficulty of correctly identifying these species away from their allopatric breeding grounds, further inhibiting knowledge of how behavior influences isolation between these species. Strong levels of postzygotic isolation are unlikely given the recent date of divergence and the long time frames required for intrinsic postzygotic isolation to accumulate in birds (Price and Bouvier 2002).

Within-Species Genetic Structure

Breeding Bicknell's Thrushes are primarily restricted to high-elevation mountaintops, but their “sky-island” distribution does not appear to limit gene flow between most populations. Our finding of limited geographic genetic structure among isolated Bicknell's Thrush populations, except for isolation by distance, supports the findings of the only previous study of geographic genetic variation by Ellison (2001), who obtained partial Control Region III sequences from 43 Bicknell's Thrush. However, we found that populations south of the St. Lawrence River and west of New Brunswick were usually similar to each other. Likewise, most populations north and east of the St. Lawrence River were not significantly different from each other but were significantly different from the southern populations. We found evidence of isolation by distance based on ND2, but all pairwise F_ST values above 0.1 compared a population north/east of the St. Lawrence to a southerly one, and our AMOVA results indicate that a significant amount of genetic variation is explained by the north vs. south groupings. Clear groupings were not supported by the SNP data, but a northern population (Nova Scotia/New Brunswick) and a southern population (northern New York) were highly differentiated from each other and most other populations. Together, these results indicate that there is slight genetic structure within Bicknell's Thrush, but the structure is not driven by their “sky-island” breeding distribution or known modern-day geographic barriers.

The northern and southern mitochondrial groupings of the Bicknell's Thrush is surprising because all populations except Nova Scotia are <200 km from their nearest neighbor and most Bicknell's Thrushes disperse less than 200 km from their natal grounds (Studds et al. 2012). However, these genetic groupings may correspond to morphological differences within Bicknell's Thrush; Wallace (1939) described 2 color morphs, with “brown-backed” individuals found only in the southern part of their range, and an “olive” morph found in the Canadian Maritimes and north of the St. Lawrence River. Wallace's sample size was small and he found many “intermediates or intergrades” (p. 234). Further study of phenotypic variation within Bicknell's Thrush should examine mandibular coloration in addition to plumage coloration, as Gray-cheeked Thrush subspecies are distinguishable by the extent of yellow on the lower mandible (FitzGerald et al. 2017, Whitaker et al. 2018). We note, however, that ancestry analyses using 5,633 SNPs did not distinguish significant groupings, indicating that regular admixture occurs between populations.

As discussed in FitzGerald et al. (2017), mitochondrial data and morphological analyses support the current taxonomic splitting of Gray-cheeked Thrushes into 2 subspecies: C. m. minimus in Newfoundland and southeastern Labrador, and C. m. aliciae from central Labrador west to Siberia, although the subspecies were not reciprocally monophyletic. Interestingly, in the current study, ancestry analyses from SNP data revealed that individuals from western Labrador were more similar to C. m. minimus than to individuals from Alaska/Siberia, but segregated at K = 3; however, the most likely number of groupings was ambiguous. This indicates that genetic differentiation between the subspecies is slight and probably only detectable at quickly evolving loci, but the rates of mutations at these nuclear SNPs are unknown. Nuclear and mitochondrial genes have different mutation rates and unique inheritance patterns (Toews and Brelsford 2012), potentially explaining why Gray-cheeked Thrush subspecies showed higher support with mtDNA vs. nuclear SNPs in this study.

Effects of LGM Refugia on Species Boundaries and Population Structure

Our analyses indicate that Bicknell's Thrush and Gray-cheeked Thrush largely resided in allopatric refugia, with the potential for up to 8.2% or 12.2% overlap predicted depending on the LGM scenario. Bicknell's Thrush refugia were smaller and more southerly than those of the Gray-cheeked Thrush (Figure 4), reminiscent of their modern-day breeding ranges. Bicknell's Thrush is currently found in regions with a high abundance of balsam fir, whereas the Gray-cheeked Thrush is associated with a broader range of habitats, often including black spruce as well as alders, willows, and other deciduous shrubs (FitzGerald 2017, Whitaker et al. 2018), and the SDMs show that balsam fir and black spruce are also important contributors for the thrush models, respectively (FitzGerald 2017). Both spruce and fir expanded northward as the ice receded, but pollen records indicate that spruce migrated more quickly, reaching areas ∼1,000 yr before fir (Davis 1983). In the modern era in which climate change is largely responsible for driving species range shifts, deciduous shrubs are colonizing the boreal–arctic ecotone even more quickly than conifers (Myers-Smith et al. 2011), coinciding with a Gray-cheeked Thrush expansion in northern Labrador (Whitaker 2017) and Alaska (Mizel et al. 2016). If Bicknell's Thrush primarily tracked habitats that included balsam fir while Gray-cheeked Thrush followed black spruce or shrubs, these 2 thrush species may have expanded northward at different times or rates, helping to maintain divergence in the past and continuing into the present day.

LGM models (Figure 4) show substantially less extensive suitable habitat compared to their modern breeding ranges (Figure 1A inset), and Tajima's D and Fu's F calculated from ND2 were significantly negative for both Bicknell's Thrush and Gray-cheeked Thrush (Table 1), supporting a recent population expansion that likely occurred as boreal habitat expanded. The MIROC LGM model revealed a single, contiguous refugium in the southern United States for Bicknell's Thrush, whereas the CCSM4 scenario projects 2 distinct refugia at the LGM, consistent with the northern and southern population structure in our DNA data. Gray-cheeked Thrushes may have resided in multiple refugia during the LGM based on our Pleistocene-projected SDMs, and we speculate that the Newfoundland Gray-cheeked Thrush may have diverged while isolated in a refugium on the Grand Banks off present-day Newfoundland. The Newfoundland population may have retained its genetic distinctiveness from the Northern subspecies in part because of a distributional gap between the subspecies, currently ∼200 km wide (FitzGerald et al. 2017). Population structure cannot be assessed in detail in the central or western parts of the Gray-cheeked Thrush breeding range without more samples, but genome-wide F_ST and STRUCTURE analyses do suggest that Siberian and Alaskan samples are genetically differentiated from Labrador samples, leaving the possibility that these far-west populations might have diverged in a Beringian refuge.

The Bicknell's Thrush and Gray-cheeked Thrush diverged in the mid- to late Pleistocene (Weir and Schluter 2004, Voelker et al. 2013, this study). Although gene divergence is expected to predate species divergence, these estimates are much older than the LGM scenarios we explored with our distribution models. Climate layers reflecting environments older than 150 ky do not currently exist (see  www.worldclim.org), precluding exploration of the species distributions at the time of divergence using SDM methods. We also caution that SDMs presume niche conservatism, which may not be applicable over hundreds of thousands of years.

Performance of Pleistocene Scenarios

Validating the locations of Pleistocene bird refugia is difficult in the absence of avian fossils, but plant macrofossil and pollen records combined with genetic data from modern populations have revealed a fairly complete picture of range dynamics in tree species since the LGM. Black spruce is composed of several distinct genetic clusters, and palynological data and macrofossils show that black spruce may have resided in 5 disjunct refugia at the LGM: Beringia, Oregon/Washington, the central United States, the mid-Atlantic Coast/Appalachia, and on the Atlantic shelf off Newfoundland (Davis 1983, Jackson et al. 1997, Jaramillo-Correa et al. 2004, Williams et al. 2004, Brubaker et al. 2005). Our Gray-cheeked Thrush SDM from the CCSM4 climate scenario (Figure 4) revealed similar refugia to the black spruce SDM for the CCSM4 scenario ( Supplemental Material Figure S3). The MIROC scenario showed a large swath of suitable habitat in the southern/central United States and the mid-Atlantic Coast, but hindcasted no suitable habitat in western North America or Newfoundland. Whereas spruce was widespread throughout North America south of the Laurentide ice sheet, the pollen and macrofossil records dating from 15–21 kbp shows that species of fir (Abies) occurred in patchy, small, and isolated regions in the south/central United States and off the mid-Atlantic Coast (Davis 1983, Jackson et al. 1997). Both LGM scenarios identified these regions as potential refugia for balsam fir, and the CCSM4 model also predicted suitable habitat off Newfoundland; Bicknell's Thrush LGM refugia (Figure 4) were similar to that for balsam fir ( Supplemental Material Figure S3). Our work stresses the need to examine different modeling scenarios and to attempt to validate models using multiple lines of evidence when available.

What Makes the Evolutionary History of the Bicknell's/Gray-Cheeked Thrush Species Complex Unique?

Many boreal bird species complexes have congruent biogeographic breaks in western North America, a pattern that is often interpreted as vicariant speciation in allopatric ice-free forest refugia (e.g., Zink 1994, Weir and Schluter 2004, van Els et al. 2012). Some studies have shown that bird populations breeding on Newfoundland are genetically distinct, but that enough gene flow occurs between island and mainland populations to prevent speciation (Colbeck et al. 2008, van Els et al. 2012, Lait and Burg 2013, Dohms et al. 2017, FitzGerald et al. 2017). The unique (among birds) eastern species-level break between Bicknell's Thrush and Gray-cheeked Thrush is similar to the break between black spruce and red spruce (P. rubens). Red spruce is found primarily in cool, well-drained soils in the Appalachians, New York, New England, southern Québec, New Brunswick, and Nova Scotia, whereas black spruce is found in abundance from Alaska across the continent to the Atlantic in diverse plant communities and soil types. Genetic evidence, however, shows that red spruce is phylogenetically nested within black spruce (Perron et al. 2000, Jaramillo-Correa and Bousquet 2003), a pattern we do not see in our thrush species complex. This suggests that either (1) cryptic species breaks in eastern North American boreal forest species may remain undiscovered or (2) another factor, such as divergent natural selection, high migratory connectivity, or interspecific competition, has played a role in the divergence of the Bicknell's/Graycheeked Thrush species complex.

If populations diverge in allopatry under similar ecological conditions, they will exhibit niche conservatism (Peterson et al. 1999, Wiens and Graham 2005, Warren et al. 2008). Because the Bicknell's Thrush and Gray-cheeked Thrush are sister taxa that breed allopatrically, yet are part of similar bird species communities in the boreal forest, we might expect that these species inhabit similar niches (Gause 1934, Svärdson 1949). Contrary to this expectation, FitzGerald (2017) found that they occupy distinct niches despite the fact that similar niche space is available to each species. In conclusion, the Bicknell's Thrush and Gray-cheeked Thrush are sister species, have a phylogeographic break that—to our knowledge—is unique among boreal bird species, and occupy divergent niches, implying that ecological selection pressures may have resulted in divergence and the maintenance of species' boundaries.