Phylogeography of a migratory songbird across its Canadian breeding range: Implications for conservation units

Abstract The objectives of this study were to describe and evaluate potential drivers of genetic structure in Canadian breeding populations of the Ovenbird, Seiurus aurocapilla. We performed genetic analyses on feather samples of individuals from six study sites using nuclear microsatellites. We also assessed species identity and population genetic structure of quill mites (Acariformes, Syringophilidae). For male Ovenbirds breeding in three study sites, we collected light‐level geolocator data to document migratory paths and identify the wintering grounds. We also generated paleohindcast projections from bioclimatic models of Ovenbird distribution to identify potential refugia during the last glacial maximum (LGM, 21,000 years before present) as a factor explaining population genetic structure. Birds breeding in the Cypress Hills (Alberta/Saskatchewan) may be considered a distinct genetic unit, but there was no evidence for genetic differentiation among any other populations. We found relatively strong migratory connectivity in both western and eastern populations, but some evidence of mixing among populations on the wintering grounds. There was also little genetic variation among syringophilid mites from the different Ovenbird populations. These results are consistent with paleohindcast distribution predictions derived from two different global climate models indicating a continuous single LGM refugium, with the possibility of two refugia. Our results suggest that Ovenbird populations breeding in boreal and hemiboreal regions are panmictic, whereas the population breeding in Cypress Hills should be considered a distinct management unit.


| INTRODUCTION
Over the past several decades, many species of songbirds breeding in North America have undergone substantial population declines (NABCI, 2016). In the case of declining species with extensive distributions and geographical variation in habitat selection and trends, many authors have argued for spatially explicit conservation units to ensure effective population-level conservation plans (Carter, Hunter, Pashley, & Rosenberg, 2000;Rosenberg & Blancher, 2005;Rushing, Ryder, Scarpignato, Saracco, & Marra, 2016). However, many different population concepts have been used to delineate conservation units for migratory songbirds (e.g., global, breeding, regional populations, and subpopulations) and clearer guidance based on genetic and demographic data is required to inform conservation efforts (e.g., Moritz, 1994;Rushing et al., 2016).
In theory, one of the most direct approaches to identify conservation units is to find individuals sharing unique phenotypic and genotypic attributes in a defined area. If these unique attributes exist, that would suggest that the population is locally adapted and regionallyshared phenotypic attributes and subtle morphological characteristics have led to the recognition of subspecies (Patten & Unit, 2002;Phillimore & Owens, 2006). It has been argued that genetic assessments of subspecies should focus on putatively adaptive markers rather than neutral markers to distinguish among phenotypes (Patten, 2015; but see Moritz, 1999;Manthey, Klicka, & Spellman, 2011).
However, genetic structure analyses across species ranges have been conducted for only a few migratory songbird species (e.g., Colbeck, Gibbs, Marra, Hobson, & Webster, 2008;Ralston & Kirchman, 2012) and use of molecular markers coding for specific phenotypic traits is often not an option for most nonmodel species (but see Van Bers et al., 2012;Ruegg et al., 2014).
Genetic structure, isolation, and drift can all result from historical processes. For example, Pleistocene glaciations are thought to have led to vicariance and subsequent speciation among northern avifauna (Johnson & Cicero, 2004;Lovette, 2005;Weir & Schluter, 2004). Although a single contiguous boreal forest glacial refugium in the southeastern United States is generally recognized (Dyke, 2005;Jackson & Overpeck, 2000;Overpeck, Webb, & Webb, 1992), wide-ranging boreal-breeding bird species could have also had multiple geographically isolated refugia across North America during the Pleistocene epoch, leading to genetic and behavioral divergence.
Such historical isolation may explain current differences in migratory paths and wintering areas. However, the importance of historical versus current barriers and contemporary demographic processes on limiting migration and dispersal pathways-and, ultimately, influencing the persistence of distinct genetic structure-is relatively unexplored (but see Ruegg et al., 2006;Milá, McCormack, Castaneda, Wayne, & Smith, 2007;Delmore, Fox, & Irwin, 2012). Thus, the degree to which populations are demographically (social interactions via natal or breeding dispersal) or genetically isolated is effectively unknown for most migratory songbirds.
For migratory species, understanding the geographical connections between different phases of the annual cycle (i.e., migratory connectivity) is essential for effective conservation (Marra, Hunter, & Perrault, 2011). When migratory connectivity is strong, wintering and breeding populations are considered as discrete units. Weak migratory connectivity implies considerable mixing among populations on the wintering grounds or large-scale dispersal among breeding areas (Webster, Marra, Haig, Bensch, & Holmes, 2002). Strong migratory connectivity makes it easier to assess the efficacy of full life-cycle management actions (e.g., Taylor & Stutchbury, 2016). With the development of tracking technologies, migratory connectivity can now be quantified using light-level geolocation, radiotelemetry, or GPS tags Stutchbury et al., 2009;Woodworth, Mitchell, Norris, Francis, & Taylor, 2015).
Even minimal movement among bird populations can result in genetic panmixia (e.g., Cortes-Rodriguez, Sturge, & Omland, 2016), further complicating the distinction of conservation units. However, cryptic genetic structure in such apparent panmixia can often be identified using parasites, mutualists, and commensals. Symbionts of birds have been shown to be useful biological markers for resolving host evolutionary history, population structure, and migratory patterns because of their shorter generation time compared to their hosts (Bruyndonckx, Biollaz, Dubey, Goudet, & Christe, 2010;Palopoli et al., 2015;Whiteman, Kimball, & Parker, 2007). Variation in symbiont assemblages or in genetic structure of particular species of symbionts can potentially provide finer resolution of host structure than the hosts themselves. However, not all symbionts are equally likely to reflect host population structure. Quill mites (Acariformes, Prostigmata, Syringophilidae) may be good indicators of host population structure, as they have no free-living stage, undergo several generations per year on one host (Kethley, 1971), and are transmitted among birds only through close physical contact (e.g., between mating individuals, parent-offspring, among nestlings). Unlike some lice and skin mites (Epidermoptidae), quill mites have never been reported as phoretic on vagile parasitic flies (e.g., Harbison & Clayton, 2011). We used intrinsic and extrinsic markers to examine historical processes that may be responsible for the contemporary distribution of Ovenbird in Canada. Specifically, we conducted genetic analyses of birds and of their syringophilid quill mites sampled at six study sites.
Results from genetic analyses were compared to those from light-level geolocators deployed on individuals from populations in eastern and western Canada. In doing so, we evaluated the hypothesis that there would be population structure in Ovenbird across its Canadian breeding range. Population structure could occur among migratory flyways, whereby individuals from the western portion of the species' Canadian breeding range would show strong migratory connectivity and be more genetically (host and symbionts) similar than those breeding in the eastern portion of the breeding range, despite the absence of physical barriers. Population structure could also occur both within and among migratory flyways as a result of finer-scale migratory connectivity (e.g., leapfrog migration; Boulet & Norris, 2006). Alternatively, there could be enough dispersal movements (breeding and natal) among Ovenbirds in Canada to prevent genetic differentiation (host and symbionts; i.e., panmictic population). These alternatives were evaluated against paleoclimatic hindcasts of potential LGM refugia. Using these multiple lines of inquiry, we aimed to evaluate evidence for Ovenbird population structure that could in turn provide guidance for efficient population management and conservation.

| Study area and data collection
In 2012, territorial male Ovenbirds were captured in the Black Brook
Full details of microsatellite development are in Appendix S1. All 15 microsatellites were scored using GeneMapper 4.0. Microchecker 2.2 (Van Oosterhout, Hutchinson, Willis, & Shipley, 2004) was used to check for errors in the data and presence of null alleles. To ensure the neutrality of markers, we tested for linkage disequilibrium and deviation from Hardy-Weinberg equilibrium using genepop 4.2 (Rousset, 2008). Unbiased expected (H E ) and observed heterozygosity (H O ) were calculated using Msa 4.05 (Dieringer & Schlötterer, 2003).
A principal coordinates analysis (PCoA) was performed in genalex 6.5 (Peakall & Smouse, 2012) using G′ ST estimates (Hedrick, 2005) to visualize population differentiation. To assess the level of population genetic structure present, we used structure 2.3 (Pritchard, Stephens, & Donnelly, 2000) with the following parameters: 10 iterations per K F I G U R E 1 Study sites where feather samples from male Ovenbirds were collected and/or geolocators deployed (a). Also indicated are the breeding and wintering ranges according to BirdLife International and Nature Serve (2012). The red line depicts the median locations during spring migration and the light red error distribution represent 95% credible interval of the migration locations (b-f). The migratory routes and wintering grounds of S. a. aurocapilla individuals fitted with geolocators are indicated with a red color ramp. The map projection is North America Albers equal area (i.e., clusters), a total of 500,000 MCMC and 50,000 burn-in, admixture model, correlated allele frequencies, and testing K 1-8. STRUCTURE analyses were performed with and without location prior information (LOCPRIOR = 0 or 1; Hubisz, Falush, Stephens, & Pritchard, 2009) to improve detection of weak population structure. We used STRUCTURE HARVESTER (Earl & vonHoldt, 2012) to determine the optimal K, using both the natural log likelihood, Ln Pr (K), and second-order rate of change of likelihood, ΔK (Evanno, Regnaut, & Goudet, 2005).
To complement assessments of population structure, we estimated the directional component of genetic divergence using DivMigrate (Sundqvist, Keenan, Zackrisson, Prodohl, & Kleinhans, 2016). To account for uneven sample sizes, we randomly subsampled from larger populations 10 times. The resulting 10 data sets comprised 22 individuals in each population, equal with the smallest sample size (Thunder Bay). DivMigrate was run with a 0.05 alpha, calculating G ST statistics for network plots and displaying connections >0.5 in strength. G ST matrix values were averaged across the 10 data sets.

| Mite identification and genetic analyses
The calamus of each feather sample was examined for syringophilid mites using a dissecting microscope (see Fig. S2.1 in Appendix S2).
In addition to the 363 feather samples used for bird genetic analyses, more samples from both LSLBO and Black Brook were examined.
Calami with mites were cut and stored in 95% ethanol. DNA extractions were performed on calami with the mites in situ. If both calami from a single host contained mites, they were placed in the same tube and extracted as a single sample. See Appendix S2 for a detailed description of the amplification of a fragment of the cytochrome oxidase subunit I (COI) gene and technique used to slide mount 384 syringophilids from 20 host birds.

| Light-level geolocator analyses
In 2012, 13 individuals were fitted with archival light-level geolocators in Black Brook and LSLBO (hereafter "geolocator," 0.65 g, Mk20S, Lotek Wireless, Newmarket, Canada) (n = 26). Each geolocator was glued to a leg-loop harness that was fitted to a bird (after Rappole & Tipton, 1991). One or three color bands (one or two bands per leg) were also applied to each individual fitted with a geolocator for visual identification in the subsequent year based on unique color combinations.
In 2013, we retrieved four units from Black Brook and one from LSLBO, for a retrieval rate of 19%. We also deployed 20 geolocators at three additional locations: (1) LSLBO, (2) Fort McMurray, and (3) Cypress Hills Interprovincial Park (a total of 60 geolocators were deployed). In 2014, we retrieved seven geolocators (two from LSLBO, two from Fort McMurray, and three from Cypress Hills), for a retrieval rate of 12% (14% across the 2 years). In total, 19% (16/86) of individuals marked in year x returned in year x + 1 (i.e., five birds could not be recaptured, loss their harness or geolocators failed to record data).
Light data collected by the geolocators were converted to geographical locations (latitude and longitude) to document vernal migration and identify wintering grounds using the Solar/Satellite Geolocation for Animal Tracking package ("SGAT"; Wotherspoon, Sumner, & Lisovski, 2013;Sumner, Wotherspoon, & Hindell, 2009) (Daly et al., 2008). We used this model to project refugia based on paleoclimate projections for the last glacial maximum (LGM).

| Paleohindcast models
As inputs to the distribution model, we sampled 10% of all 4-km pixels across North America (n = 139,890) and intersected them with the range map to assign the presence/absence locations (prevalence = 0.256) as the dependent variable and with derived bioclimatic variables as independent predictors. We calculated seven bioclimatic variables chosen based on relevance to vegetation distributions (Hogg & Bernier, 2005), avoidance of extreme collinearity (Dormann et al., 2013), and a preference for seasonal over annual variables when they showed high correlations (Stralberg et al., 2015).
These variables included extreme minimum temperature, chilling degree days, growing degree days, seasonal temperature difference, mean summer precipitation, climate moisture index, and summer climate moisture index (available as a data supplement to Stralberg et al., 2015).
We used the "dismo" (Hijmans, Phillips, Leathwick, & Elith, 2011) and "gbm" (Ridgeway, 2012) packages for R (R Core Team, 2014) to build a distribution model for the species, specifying a binomial distribution ("bernoulli" family). We used a stepwise procedure and 10-fold crossvalidation to identify the optimal number of trees needed to maximize the mean deviance explained as recommended by Elith (2) the Geophysical Fluid Dynamics Laboratory (GFDL) model, from the National Oceanic and Atmospheric Administration. The monthly anomalies with respect to baseline climate variables (Roberts & Hamann, 2015; In press) were used as inputs to the BRT models to develop millennial-scale hindcasts using the "raster" package for R (Hijmans & van Etten, 2012). Core habitats (both current range and LGM refugia) were defined as areas of predicted occupancy greater than 0.256 within the area covered by the training dataset, corresponding to the prevalence rate.
Finally, to compare climatic niches among subspecies, we repeated the above-described process for mapped subspecies S. a. aurocapilla, S. a. furvior, and S. a. cinereus (adapted from American Ornithological Union, 1957). We also ran a principal component analysis on the three subspecies, using the "prcomp" function for R based on centered and scaled versions of the seven bioclimatic variables used to develop distribution models.

| Bird genetics
Of the 363 feather samples collected, 339 were successfully genotyped at 12 or more loci. The presence of null alleles and linkage disequilibrium was not detected. Eight loci remained after the removal of seven loci deviating from HWE (see Table S3.2 in Appendix S3).
The total number of alleles per locus ranged from 16 (Dpu16) to 31 (Dpu01). Population H E was lowest for Cypress Hills and Thunder Bay (0.76) and highest for Gatineau and Barachois Pond (0.80; Table 1).  S1.3 in Appendix S1). When Cypress Hills was removed from the STRUCTURE analysis, we found no evidence for further substructure despite retaining sampling location as a prior. DivMigrate suggested a similar pattern in which gene flow was present among all populations except Cypress Hills (Figure 3). Individual network plots for the 10 subsampled data sets are in Appendix S1.  F I G U R E 2 structure assignment of individuals to each genetic cluster for K = 2 using eight microsatellites. Population structure was weak when sampling locations were not included as a prior in the analysis (a) and only detected when sampling locations were included (b) distinct clustering in host populations and very little variation among specimens. Since the mean genetic distance among the Cypress Hills specimens is larger than the mean distances between the Cypress Hills mites and those from other areas (Table S2.1 in Appendix S2), the Cypress Hills mites cannot be considered a distinct taxon (Čandek & Kuntner, 2015).

| Light-level geolocator analysis
Geolocators recovered from New Brunswick (Black Brook) suggest that individuals from this population overwintered in Cuba, Hispaniola, or the Bahamas (Figure 1b

| Paleohindcast models
Predictive accuracy of the BRT distribution model was high: 94% as measured by the 10-fold crossvalidation R 2 . Spatial predictions from paleohindcasting based on both climate models suggest that the Ovenbird likely had a single contiguous refugium in the southeastern United States, overlapping with the southern limit of its current range, and possibly extending more patchily through the southwestern interior United States ( Figure 5). According to the GFDL model, some suitable climate space was also available in California.
A principal components analysis revealed high overlap in climatic niche space for the three subspecies, with S. a. furvior com-

| DISCUSSION
Our results provide little evidence for genetic differentiation among Ovenbird populations in Canada, with the exception of the southwestern Cypress Hills population (S. a. cinereus). There was also limited genetic variation among symbionts (syringophilid mites) collected from all populations. These results are consistent with paleohindcast predictions from two different GCM, suggesting that a fairly contiguous glacial refugium was available in the southeastern United States 21,000 years BP, with the possibility of a separate western refugium (see also Stralberg et al., In press). We did, however, find support for differences in migratory patterns between eastern and western populations and those from northern versus southern Alberta, with evidence of relatively strong migratory connectivity within each of these three populations. The distinct genetic structure for S. a. cinereus, which could be maintained by the strong migratory connectivity exhibited by this population, provides some support for the current range delineation for this subspecies. However, no such evidence was found in support of S. a. furvior. Taken together, these results indicate that to be relevant, management units for migratory birds must be defined on the basis of detailed genetic data, but also information on migratory connectivity and glacial refugia.  (Cobben et al., 2011;Dyke, 2005;Hewitt, 2000). However, this level of genetic admixture likely contributes to fuzzy phenotypical and genetical boundaries among subspecies. Genetic information from a larger number of populations would be required for a more detailed assessment of the phylogeography of our focal species. It is also important to consider that the weak genetic structure we observed may reflect the resolution of markers used (e.g., Kimura et al., 2002;Ruegg et al., 2014).

| Lack of population structure in symbionts
Syringophilids were found on at least one bird from each population, and all the specimens belonged to the same species (B. seiuri), resulting in no geographical variation in mite assemblages. There was also no genetic variation in syringophilids that was correlated with host localities. This may reflect host populations mixing on the wintering grounds or the more extensive breeding and natal dispersal movements generally observed in female songbirds compared to males (Greenwood, 1980;Greenwood & Harvey, 1982). It is also possible that there was genetic variation associated with host localities, but the marker used (COI) might not have reflected this. Other symbiont taxa might provide a higher resolution in host population structure.
For example, members of the vane-dwelling feather mite genera Trouessartia and Proctophyllodes were also collected from some rectrices (H. Proctor & A. Grossi, unpublished data), but at numbers too low to use as markers in this study. These mites may be more abundant on other feather types, and future studies could include their taxonomic and genetic analysis.

| Migratory connectivity and population structure
With the exception of one individual breeding in Alberta that overwintered in either Florida or Cuba, migratory connectivity was strong, which is consistent with results from . Few studies have documented weak migratory connectivity or population mixing in songbird species (e.g., Ruegg et al., 2014;Hobson & Kardynal, 2015; but see Finch, Butler, Franco, & Cresswell, 2017).
Future studies should investigate the demographic implications of the relatively low mixing (9% of individuals; 1/11) reported in this study.
Although geolocator data provide important information on species distribution on breeding and wintering grounds and their migratory routes, units have only been deployed on adult males (i.e., after hatchyear) owing to higher site fidelity than females and hatch-year birds.
The implications of natal dispersal movements on the strength of migratory connectivity reported for songbirds (adults) remain unknown (Cresswell, 2014) and can only be quantified using other approaches such as stable isotopes, genetic analyses, or GPS tracking (e.g., Clegg et al., 2003;Kelly, Ruegg, & Smith, 2005).
Some studies have reported leapfrog migration in songbirds, where individuals breeding in the northern portion of the breeding range winter further south than southern-breeding individuals (Bell, 1997;Boulet & Norris, 2006;Kelly, Atudorei, Sharp, & Finch, 2002). This migration system is consistent with the pattern observed in Ovenbirds Birds from Cypress Hills may also use the Pacific flyway, as individuals are being observed somewhat regularly in California during migration (Porneluzi et al., 2011). Different migration strategies among populations could lead to morphological differentiation (Winkler & Leisler, 1992; see also Desrochers, 2010 song characteristics should also be included in a more rigorous analysis of phenotypes across the breeding range (Patten, 2015).

| Conservation implications
Population targets are being set via species-level conservation plans such as Recovery Strategies and Action Plans for species at risk This hierarchical organization of individuals could be considered a simple concept applicable to a broader range of taxa, but such spatiallyexplicit information is only available for a few migratory songbird species. Conservation efforts should focus on maintaining genetic diversity or evolutionarily significant units (e.g., Moritz, 1994Moritz, , 1999Palsboll, Bérubé, & Allendorf, 2006) and population dynamics through reproduction, survival, and functional/migratory connectivity (e.g., Rushing et al., 2016). Our results suggest that subspecies defined by arbitrary criteria do not provide a strong model for setting population management units. For example, conservation of Ovenbirds in Newfoundland (S. a. furvior) should be based on a larger area given the extensive gene flow that occurs with S. a. aurocapilla. In contrast, a subspecies-based conservation strategy would be appropriate for S. a. cinereus in southern Alberta and Saskatchewan. Subspecies have been identified based on similar criteria for species of conservation concern (e.g., Rusty Blackbird, Euphagus carolinus nigrans, American Ornithological Union, 1957; and White-winged Crossbill, Loxia curvirostra percna, Dickerman, 1987), and genetic analyses are required to inform population management plans. We encourage future studies such as this that employ multiple lines of evidence to identify different levels of population structure for conservation purposes.

ACKNOWLEDGMENTS
This study was supported by discovery grants from the Natural