Genetic structure and biogeographic history of the Bicknell's Thrush/Gray-cheeked Thrush species complex

ABSTRACT We examined species limits, admixture, and genetic structure among populations in the Bicknell's Thrush (Catharus bicknelli)–Gray-cheeked Thrush (C. minimus) species complex to establish the geographic and temporal context of speciation in this group, which is a model system in ecology and a high conservation priority. We obtained mitochondrial ND2 sequences from 186 Bicknell's Thrushes, 77 Gray-cheeked Thrushes, and 55 individuals of their closest relative, the Veery (C. fuscescens), and genotyped a subset of individuals (n = 72) at 5,633 anonymous single nucleotide polymorphic (SNP) loci. Between-species sequence divergence was an order of magnitude greater than divergence within each species, divergence was dated to the late Pleistocene (420 kbp) based on Bayesian coalescence estimation, and a coalescent model (IMa) revealed almost no gene flow between species based on ND2. SNP data were consistent with mitochondrial results and revealed low levels of admixture among species (3 of 37 Bicknell's Thrushes, no Gray-cheeked Thrushes, and no Veeries were >2% admixed). Species distribution models projected to the Last Glacial Maximum suggest that Bicknell's Thrush and Gray-cheeked Thrush resided in primarily allopatric refugia in the late Pleistocene, consistent with the genetic data that support reproductive isolation over an extended period of time. Our genetic data suggest that both species underwent demographic expansions, possibly as they expanded out of Pleistocene refugia into their current ranges. We conclude that Bicknell's Thrush and Gray-cheeked Thrush are 2 distinct species-level lineages despite low levels of Gray-cheeked Thrush introgression in Bicknell's Thrushes, and divergence has been maintained by a long history of allopatry in subtly different habitats. Their unique phylogeography among boreal forest birds indicates that either cryptic species breaks in eastern North America are still undiscovered, or another factor, such as divergent natural selection, high migratory connectivity, or interspecific competition, played a role in their divergence.


INTRODUCTION
Most speciation events are suspected to occur in allopatry. Reproductive isolation can arise when single mechanisms of divergence, such as natural selection to distinct ecological conditions, sexual selection, or genetic drift, act in conjunction with geographic isolation (Coyne andOrr 2004, Price 2008). Even for populations occupying large contiguous landmasses, allopatric speciation is predominant, as supported by evidence of niche conservatism between sister species (Peterson 2011). Similarly, congruent biogeographic breaks in co-distributed species and species complexes imply that multiple lineages were split due to a single event or barrier, although that barrier may no longer be present Kenagy 2001, Coyne andOrr 2004). For example, during the Pleistocene, glacial ice sheets expanded southward and inundated northern habitats above ~42°N during cool periods and retreated northward when temperatures warmed (Hewitt 2004). Although divergence in ice-free allopatric refugia has been questioned on the basis of incongruent timing Zink 1997, Arbogast andSlowinski 1998), many studies have found a detectable genetic signature of isolated evolution among populations that were sundered by Pleistocene ice sheets (e.g., Hewitt 2004, Weir and Schluter 2004, Lait and Burg 2013, Dohms et al. 2017. The modern-day boreal forest biome of North America was entirely ice-covered during the Last Glacial Maximum (LGM). Species and species complexes that breed today in the boreal forest were pushed out of the area into smaller, disjunct refugia separated by uninhabitable regions during the LGM, and therefore show stronger signatures of isolation in glacial refugia compared to temperate and tropical species (Weir and Schluter 2004). Surveys of genetic population structure in boreal forest groups have found broadly congruent patterns of divergence between a single "boreal" clade in northern and eastern North America and one or more clades west of the Rocky Mountains, with most of the genetic structure (and morphological differences, if present) found in the western portions of species' ranges (Zink 1994, Weir and Schluter 2004, van Els et al. 2012, Dohms et al. 2017). This pattern is exemplified by the Canada Jay (Perisoreus canadensis) species complex, which comprises a genetically diverse "boreal" clade stretching from eastern North America to Alaska, and multiple, lower-latitude western clades; these clades may represent separate species based on genetic distances, although they are not currently recognized as such taxonomically. Species distribution models (SDMs) of Canada Jays at the LGM identify 4+ western refugia and 1-3 eastern refugia (van Els et al. 2012, Dohms et al. 2017. Boreal species like Blackpoll Warbler (Setophaga striata) and Boreal Chickadee (Poecile hudsonicus) are not found southwest of the Rockies and do not have deep biogeographic breaks across their range. These species have similar geographic structures, with imperfect clustering of western, eastern, and Newfoundland birds and high levels of gene flow between clusters Kirchman 2012, Lait andBurg 2013), and their LGMprojected SDMs indicate possible refugia in Beringia, the southeastern United States, and Newfoundland. Pollen and fossil records of boreal tree species suggest that a large, relatively contiguous boreal-type forest refugium existed in the southeastern United States, multiple, smaller refugia persisted in Beringia and western North America, and a small sub-boreal refugium could be found on the then exposed continental shelf off the coast of Newfoundland during the LGM (Davis 1983, Jackson et al. 1997, Jaramillo-Correa et al. 2004, Williams et al. 2004, Brubaker et al. 2005, largely corroborating the patterns found in birds. Together, these studies indicate that (1) boreal bird species and species complexes diversified during the Pleistocene; (2) interspecific-level breaks are found in western North America, whereas there is little geographic genetic structure across the broad swath of boreal forest in central and eastern North America excepting Newfoundland; and (3) SDMs projected to the LGM support within-species patterns of genetic structure, suggesting that modern genetic breaks are the result of geographic isolation during the Pleistocene.
Among North American boreal forest birds, only one species complex has a species-level break in eastern North America: the Bicknell's Thrush (Catharus bicknelli Ridgeway) and Gray-cheeked Thrush (C. minimus Lafresnaye) ( Figure 1A inset). From its earliest descriptions, the Bicknell's Thrush was considered a subspecies of Gray-cheeked Thrush (Ridgway 1881, Bicknell 1882 American Ornithological Society A. M. FitzGerald, J. Weir, J. Ralston, et al. Phylogeography of Catharus bicknelli/C. minimus complex Brewster 1883, Wallace 1939, and was not recognized as a full species by taxonomic authorities (e.g., Monroe et al. 1995) until Ouellet (1993) showed that it is distinct from C. minimus in vocalizations, morphology, breeding habitat, and wintering range. However, Marshall (2001) argued that, because songs were very similar and bodysize and plumage variation were clinal, the Newfoundland subspecies of Gray-cheeked Thrush (C. m. minimus) was a geographical and phenotypic intermediate between the  (Dilger 1956, Able and Noon 1976, Noon 1981) and modes of speciation (Winker 2010, FitzGerald 2017, has been the focus of several molecular phylogenetic studies aimed, in particular, at understanding the evolution of migration (Ruegg and Smith 2002, Outlaw et al. 2003, Winker and Pruett 2006, Voelker et al. 2013. Previous studies have shown that, together with the Veery (C. fuscescens), Gray-cheeked Thrush and Bicknell's Thrush comprise the most recently diverged lineage within the genus, but these studies have inferred different sister species relationships within this clade depending on the number and type of genetic loci (i.e. mitochondrial vs. autosomal) sequenced (see Ellison 2001, Outlaw et al. 2003, Klicka et al. 2005, Winker and Pruett 2006. The most recent and comprehensive phylogenetic analysis based on sequences of 6 autosomal, 2 Z-chromosomal, and 2 mitochondrial loci from one individual per species found high support for Gray-cheeked Thrush and Bicknell's Thrush as sister species with the Veery as sister to this clade, and estimated that this group split from its sister clade ~800,000 yr before present (ybp) (Voelker et al. 2013).
The only genetic studies to sample multiple individuals (i.e. more than 2) in the Bicknell's Thrush/Gray-cheeked Thrush species complex were Ellison (2001) andFitzGerald et al. (2017), the latter of which assessed the distinctiveness of Bicknell's Thrush and the 2 Gray-cheeked Thrush subspecies on the basis of mitochondrial and nuclear DNA sequences and morphology. FitzGerald et al. (2017) were unable to clearly resolve major genetic breaks from 2 nuclear loci but found that mitochondrial divergence between Gray-cheeked Thrush subspecies was an order of magnitude smaller than the divergence between Graycheeked and Bicknell's thrush, supporting a species-level break. However, they also reported the first putative hybrid between the Gray-cheeked Thrush and Bicknell's Thrush, evidence of occasional introgression. The morphological and genetic results supported designation of Gray-cheeked Thrushes from the Newfoundland archipelago and Harbour Island, Nova Scotia, as C. m. minimus, and birds from western Labrador and farther west as C. m. aliciae, with southern Labrador comprising a zone of mixing between the subspecies .
Lack of concordance between mitochondrial and nuclear gene trees in the recently diverged Veery/Gray-cheeked/ Bicknell's thrush lineage may stem from recent hybridization or incomplete lineage sorting. Examining many individuals per species at many loci sampled throughout the genome can provide stronger statistical power for resolving phylogenies and enable rigorous tests of genetic admixture in closely related species complexes (e.g., Harvey and Brumfield 2015, Mason and Taylor 2015, Winger et al. 2015, Jeffries et al. 2016). Here we implement genotyping-bysequencing, a restriction-site associated DNA sequencing method that results in thousands of polymorphic loci (Elshire et al. 2011), along with Sanger sequencing of mitochondrial DNA from the largest sample size of thrushes to date. We combine these genetic analyses with Pleistocene distribution models to investigate the biogeographic history of the Bicknell's/Gray-cheeked Thrush species complex.
Our analysis of genetic structure and biogeographic history provides a new basis for taxonomic treatment of the group, but also provides important evolutionary context to the ecological study of these threatened and declining taxa. Currently, the Bicknell's Thrush is considered the "Northeast's most vulnerable songbird" (Rimmer and McFarland 2013) because of its small and naturally fragmented breeding range, its small and vulnerable wintering range (Lambert et al. 2005, Rodenhouse et al. 2008, COSEWIC 2009Matteson 2012, Environment Canada 2014, and documentation of significant declines in the last 2 decades (Ralston et al. 2015). Because of this conservation concern, Bicknell's Thrush is the primary focus of the Mountain Birdwatch citizen-science monitoring effort (vtecostudies.org; Rimmer and McFarland 2013), which has resulted in a precise estimate of its global population size of <120,000 birds (71,318 [95% CI: 56,(80)(81)(82)(83)(84)(85)(86)(87)(88)(89)748] in the United States and ~45,000 in Canada; Hill and Lloyd 2017). Breeding Bird Survey (BBS) data indicate that the Newfoundland Gray-cheeked Thrush (C. m. minimus) has undergone a ~95% decline over the past 30 yr (SSAC 2010, Whitaker et al. 2015FitzGerald et al. 2017). As a result, this population has been listed as threatened under the Newfoundland and Labrador Endangered Species Act and is awaiting assessment for possible listing under Canada's federal Species at Risk Act (COSEWIC 2013). Further, demographic and BBS data indicate that Northern Graycheeked Thrushes (C. m. aliciae) in Alaska and possibly northwestern Canada are also declining, while monitoring data are virtually nonexistent for other portions of the species' range (DeSante et al. 2015, Sauer et al. 2017, Whitaker et al. 2018). Our goals were to (1) evaluate the reciprocal monophyly, levels of genetic admixture, and timing of divergence between the Gray-cheeked Thrush and Bicknell's Thrush; (2) examine patterns of population structure within each species; and (3) reveal the biogeographic context of within-species and among-species divergences at the LGM.

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  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). Highquality 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 maximumlikelihood 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).

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 medianjoining 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 andBouckaert 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 PartitionFinder 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 Graycheeked 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 ) and the maximum clade credibility tree was determined with TreeAnnotator 2.1.2 (http://beast.bio.ed.ac.uk/treeannotator) using a 25% burn-in. Node support and mean age estimates were tallied using FigTree 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 1 = θ /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 distancebased 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); 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 ) was used to assess convergence. DensiTree  was used to export trees; the Bicknell's Thrush/Graycheeked 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 lowelevation 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 Graycheeked Thrush ( Figure 1A, Supplemental Material Table  S1). Bicknell's Thrush samples from Cộte-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 A. M. FitzGerald, J. Weir, J. Ralston, et al. Phylogeography of Catharus bicknelli/C. minimus complex 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 Graycheeked 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) 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 andRotenberry 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 LGMprojected habitat suitability layers in the thrush LGMprojected 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(Warren et al. , 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 1 and q 2 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 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) 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). 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 at K = 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/ 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. Betweenspecies 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.   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. 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 ; 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; 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  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. 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 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 2 = 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. 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 2019 American Ornithological Society Phylogeography of Catharus bicknelli/C. minimus complex A. M. FitzGerald, J. Weir, J. Ralston, et al. 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 , 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 Graycheeked 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 Graycheeked Thrush (Figure 4), reminiscent of their modernday 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) 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. 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 andBousquet 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 Graycheeked 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.

ACKNOWLEDGMENTS
We wish to thank the following individuals for their help in the field obtaining blood samples: