**Prioritizing Areas for Land Conservation and Forest Management Planning for the Threatened Canada Warbler (***Cardellina canadensis***) in the Atlantic Northern Forest of Canada**

#### **Alana R. Westwood 1,\*, J. Daniel Lambert 2, Leonard R. Reitsma <sup>3</sup> and Diana Stralberg <sup>1</sup>**


Received: 14 December 2019; Accepted: 30 January 2020; Published: 4 February 2020

**Abstract:** Populations of Canada Warbler (*Cardellina canadensis*) are declining in Canada's Atlantic Northern Forest. Land conservancies and government agencies are interested in identifying areas to protect populations, while some timber companies wish to manage forests to minimize impacts on Canada Warbler and potentially create future habitat. We developed seven conservation planning scenarios using Zonation software to prioritize candidate areas for permanent land conservation (4 scenarios) or responsible forest management (minimizing species removal during forest harvesting while promoting colonization of regenerated forest; 3 scenarios). Factors used to prioritize areas included Canada Warbler population density, connectivity to protected areas, future climate suitability, anthropogenic disturbance, and recent Canada Warbler observations. We analyzed each scenario for three estimates of natal dispersal distance (5, 10, and 50 km). We found that scenarios assuming large dispersal distances prioritized a few large hotspots, while low dispersal distance scenarios prioritized smaller, broadly distributed areas. For all scenarios, efficiency (proportion of current Canada Warbler population retained per unit area) declined with higher dispersal distance estimates and inclusion of climate change effects in the scenario. Using low dispersal distance scenarios in decision-making offers a more conservative approach to maintaining this species at risk. Given the differences among the scenarios, we encourage conservation planners to evaluate the reliability of dispersal estimates, the influence of habitat connectivity, and future climate suitability when prioritizing areas for conservation.

**Keywords:** bird distribution and abundance; boreal birds; Canada Warbler; *Cardellina canadensis*; Zonation; reserve design

#### **1. Introduction**

Canada Warbler (*Cardellina canadensis*) is a Neotropical migratory songbird that breeds in forests of the eastern U.S. and across Canadian forests from Nova Scotia to the Yukon [1]. Due to ongoing population declines (71% decline reported from 1970–2010 [2]), it is listed as Threatened in Canada [3] and a Species of Greatest Conservation Need in nearly every U.S. state where it breeds (e.g., [4]). The Canada Warbler International Conservation Initiative (CWICI) has urged research and conservation actions to help reverse the population decline, including identification of suitable habitat and development of best management practices for the breeding grounds [5].

In the Atlantic Northern Forest (Bird Conservation Region (BCR) 14), Canada Warbler population declines from 1970–2017 have been much steeper than those observed in other BCRs by the Canadian Breeding Bird Survey [6]. Canada Warblers predominantly use wet forests [7–10] and post-harvest deciduous and mixed-wood forests approximately 10–30 years of age [11,12]. Management guidelines recently produced for BCR 14 [13] (see Supplementary Materials) describe two different approaches to promote recovery of this species: permanent land conservation (hereafter 'land conservation', or 'LC') and responsible forest management (hereafter 'forest management', or 'FM'). These approaches were developed in tandem with, and designed for, two different user communities with distinct needs: (1) land conservancies and government agencies with a mandate to protect species and habitat in situ; and (2) forest industry staff and government agencies with a mandate to engage in sustainable development of forest resources while minimizing impacts to migratory birds and species at risk. LC scenarios (summarized in Westwood et al. [13]) (see Supplementary Materials) are intended to support conservation activities such as in-situ permanent protection of forested areas supporting Canada Warbler populations, forest blocks with low edge-to-interior ratios, and suitable habitat patches connected by forested corridors to other potential breeding sites to facilitate dispersal. FM scenarios are intended to support responsible forest management activities such as providing a continuous current supply of breeding habitat on the landscape, avoiding harvesting and road building in population centers and forested wetlands, and locating areas to implement silvicultural systems (such as shelterwood cutting) most likely to produce desired conditions for breeding habitat 12–20 years post-harvest, among other strategies [13] (see Supplementary Materials).

Forest fragmentation has been postulated as a cause of population declines for this species [14,15], so habitat connectivity on the breeding grounds is likely to play an important role in population recovery, especially given the high degree of conspecific attraction that has been documented [16]. Climate change is another important consideration for Canada Warbler, as its range is expected to retract under warming conditions [17]. Conservation planning software makes it possible to account for connectivity, climate change, and other factors when prioritizing areas for land conservation and forest management [18–20].

Most analyses using conservation planning software are designed to support the planning or evaluation of protected area networks (e.g., [21–23]). Other objectives include locating zones of interest for further study [24], evaluating trade-offs between land uses [25,26], or estimating the value of ecosystem services [27]. Conservation planning algorithms can be informed by estimates of dispersal, which is an important factor in wildlife habitat selection and response to climate change (i.e., the ability to move from current to future suitable habitats). For the Canada Warbler, published estimates of natal dispersal (distance from an individual's birth site to their breeding site) are not available. Estimates of natal dispersal for other passerines are frequently based on small sample sizes with a great deal of uncertainty, and estimates for species of similar sizes to the Canada Warbler range widely [28–31]. Conservation planning exercises are typically based on single dispersal estimates, and the implications of uncertainty in these estimates are unknown. Furthermore, conservation planning algorithms often rely on species distribution models which predict habitat suitability, species occurrence, or population density across a landscape. However, these models are inherently more uncertain in under-sampled areas, and this uncertainty is also not always accounted for during conservation prioritization exercises.

Conservation prioritization exercises are most useful when they directly inform resource allocation decisions across landscapes [32]. Therefore, we consulted a variety of government, conservation, and forest industry stakeholders to develop regional scenarios for prioritization [33] (see Supplementary Materials). We then used the program Zonation [34] to evaluate priority areas for land conservation and forest management to support Canada Warbler populations in the Canadian portion of BCR 14. One set of scenarios was designed to prioritize areas for long-term land conservation (including future climate change and connectivity to protected areas) and another set of scenarios prioritized areas for forest management on tenured land managed by timber companies.

We evaluated three estimates of natal dispersal distance for each set of scenarios. We evaluated the impact of these factors on spatial outcomes and identified areas that were consistently prioritized. We also compared conservation efficiency (proportion of the species' current estimated population protected per unit area) across scenario types and dispersal distance estimates. The resulting rankings of the landscape in BCR14 are intended to support decisions for maintaining and managing Canada Warbler habitat in this region by the two different user groups already described.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The North American Bird Conservation Initiative (NABCI) defines BCRs as ecologically distinct regions that share similar bird communities, habitats, and resource management issues [35]. We selected the scale of the study as the Canadian portion of the Atlantic Northern Forest (BCR 14) for three reasons. First, BCRs are used by Canadian government agencies as planning and management units for other bird species at risk (e.g., [24,36]). Secondly, the Canada Warbler shows evidence of differential habitat selection across BCRs [37], and we wished to limit the modelling to a population with shared habitat requirements. Finally, limiting the scale of the study to the Canadian portion of the BCR allowed us to use interoperable data and take advantage of existing Canadian bird conservation and forestry networks comprised of government agencies, industry, and NGOs.

The Canadian portion of the Atlantic Northern Forest encompasses 20.4 million ha, including Nova Scotia, Prince Edward Island, New Brunswick, and the Gaspé Peninsula of Québec ([38]; Figure 1). This forest is within the Appalachian-Acadian Ecoregion, which also includes much of Maine and Vermont, and parts of New York and New Hampshire. Human activities in this area have influenced 98.2% of the land area [39].

**Figure 1.** The Canadian portion of the Atlantic Northern Forest (Bird Conservation Region 14) showing protected areas (green) and areas under lease for timber activity (white).

The Atlantic Northern Forest is at the interface of two biomes, and includes tree species from the northern boreal forest (e.g., Black Spruce *Picea mariana*, Eastern White Pine *Pinus strobus*, Balsam Fir *Abies balsamea*) and southern temperate deciduous forest (e.g., Red Spruce *Picea rubens*, Red Maple *Acer rubrum*, Sugar Maple *Acer saccharum*). Land ownership is predominantly private, with the remainder under stewardship of the provinces. Known areas that are permanently protected from development in this region comprise 1.4 million ha (7% of the land area), including national and provincial parks, wilderness areas, and private protected landholdings [40]. Timberlands leased or owned by industrial forestry companies cover 6.7 million ha (33% of the land area). We did not include privately owned woodlots or timberlands in our analyses as we were unable to obtain a comprehensive land use map for such properties.

#### *2.2. Scenario Development*

Recent habitat guidelines for Canada Warbler in the Atlantic Northern Forest identified two types of stewardship approaches that could benefit this species: land conservation and forest management [13] (see Supplementary Materials). Based on solicited input from land trust representatives, agency and industrial forest managers, and scientists familiar with Canada Warbler ecology, we developed four scenarios for land conservation and three for forest management (Table 1; see [33] (see Supplementary Materials) for a description of stakeholder engagement and [13] (see Supplementary Materials) for specific details on recommended responsible forest management practices for maintaining Canada Warbler populations on the landscape).


**Table 1.** Spatial prioritization scenarios developed to support land conservation and forest management planning for Canada Warbler conservation in the Atlantic Northern Forest.

#### *2.3. Zonation Spatial Conservation Prioritization*

#### 2.3.1. Data layers and Pre-Processing

To develop the seven scenarios, we acquired input data in three categories: avian, landcover, and administrative boundaries (Table 2).



Our avian datasets were obtained from the Boreal Avian Modelling (BAM) Project [41–43]. The BAM project holds the largest dataset of boreal and hemiboreal observations of birds in North America and accounts for heterogeneity in survey protocols by correcting abundance estimates for detectability [42]. To determine suspected breeding locations, we divided point occurrences of Canada Warblers into two groups (*CAWA presence 2005–2009* and *CAWA presence 2010–2015*) to capture areas showing persistent observations of Canada Warblers over time. Including recent presence information allows conservation planners to locate areas of high likelihood of extant Canada Warblers to consider for protection, and forest managers to avoid or operations in areas where they may harm, kill, or harass Canada Warbler or their nests or eggs (which are illegal activities under Canada's *Species At Risk Act, S.C. 2002*). To predict *Population density* in a 1 km grid across the study area, we used a national-scale species distribution model for Canada Warbler based on landcover and disturbance data [44]. The standard deviation of the population density estimates across multiple data subsets was used to represent *Model uncertainty*. As Zonation uses an additive approach to discount cell values by uncertainty, we chose a measure of uncertainty in the same units as the population density estimates (standard deviation), rather than using a standardized measure such as coefficient of variation. We note that no data were available on productivity or occupancy for the Canada Warbler in this region, nor fine-scale data to describe habitat characteristics (e.g., Lidar).

To account for projected climate-induced shifts in abundance, we used 4 km population density predictions from models based on climate, land-use, and topography covariates as compared to the baseline population density (detailed description in [45]). *CAWA climate baseline* included predicted population density from 1961–1990. Future projections were for the 2041–2070 time period (*CAWA climate 2050*) based upon a high-end, business-as-usual emissions scenario (A2), averaged over an ensemble of four global climate models from the Coupled Model Intercomparison Project (CMIP3) dataset [46]. Although these future projections do not incorporate anticipated lags in vegetation responses to climate change [17], we considered them a reasonable representation of long-term future habitat suitability for this species, which in this region is projected to experience relatively moderate shifts in density, as opposed to wholesale range shifts [47].

Administrative boundary layers included political and administrative borders (*Provinces*/ *states* and *BCR14*), known forestry tenures (*Working lands*: lands owned or leased by forestry companies for the purpose of industrial development; [47]), and *Protected areas* meeting any of the International Union for the Conservation of Nature protected areas classifications I-IV (lands protected for long-term biological values, [40]).

Landcover layers included *Water bodies* derived from 1:10,000 aerial imagery and layers of *Wet-poor habitat* and *Upland habitat* (derived from the Northeastern Habitat Types Mapping Initiative; [48]). All spatial data were processed using ArcGIS 10.2.2 [49]. All input and output data layers were rasters in geotiff format with a cell size of 1000 × 1000 m, projected in Canada Lambert Conformal Conic.

#### 2.3.2. Prioritization Analysis

We ran all analyses using Zonation 4.0 [18] with a mask of the land boundaries of the Atlantic Northern Forest (an inverse of the layer *Water bodies*) applied to eliminate oceans and lakes. Detailed run settings and input files are available at https://github.com/borealbirds/cawa-bcr-14. Zonation identifies areas with high concentrations of features, which are the items that are desirable to prioritize for the end use (e.g., population density, land ownership, etc.). Functions are used to apply rules to determine how the features interact and how connectivity between the features are prioritized or penalized [34].

Connectivity functions in Zonation rely on estimates of dispersal distance to determine whether features or populations are 'connected.' Due to scarce species-specific dispersal data, Carroll et al. [50] completed a prioritization analysis using an estimated dispersal distance (not specific to natal or breeding) of 10 km for landbird species. For Canada Warbler, there are no published natal dispersal estimates and only one known direct observation (10 km; L. Reitsma, unpublished data). Because natal dispersal distances of landbirds are correlated with wing length and body mass [29,30], the Canada Warbler's average size (wing span = 20–22 cm, mass = 9.5–12.5 g; Reitsma et al. 2010) suggests a median dispersal distance of 50 km [30]. However, estimates for similar-sized species vary widely, from 0.5 km to 40 km [28–31]. Betts et al. [51] measured maximum breeding dispersal distances of 1–3 km for two species of forest-dwelling warblers whose ranges overlap that of Canada Warbler (the Black-throated Blue Warbler, *Setophaga caerulescens*, and the Blackburnian Warbler, *Setophaga fusca*). To account for the uncertainty regarding dispersal estimates for this species and capture the variation in dispersal estimates for similar species, we evaluated three different natal dispersal estimates for each scenario: low dispersal distance (LDD, 5 km), medium dispersal distance (MDD, 10 km), and high dispersal distance (HDD, 50 km).

Zonation uses three primary algorithms to prioritize raster cells for selection: core area zonation (CAZ), additive benefit function, and targets [18]. In this case we did not have an a priori conservation target, which is one reason we chose Zonation over the similar software Marxan [19], which requires proportional or area targets as inputs. We chose CAZ because it ensures that the most valuable areas for each feature ("core areas") are prioritized, rather than allowing trade-offs between features, as is the case with additive benefit functions. In each iteration of the algorithm, CAZ chooses the cell with the lowest retention value as specified by the features and connectivity functions, discards it, and then recalculates the value of all remaining cells. In this way, each cell is ranked in order of its priority for selection, with the first cells selected being lowest priority and the last cells being highest priority [18].

For scenario LC1, we input the *Population density* feature (weight 1.0). In Zonation, features are generally given a standard weight of 1.0 unless they are to be discounted due to uncertainty (such as future climate projections) or increased in value due to particular management objectives [18]. We then applied the 'Distribution Smoothing' function to aggregate areas with cells of high population density connected by dispersal ability of Canada Warbler (dispersal kernel specified by α [18]; LDD = 5 km, <sup>α</sup> <sup>=</sup> <sup>4</sup> <sup>×</sup> <sup>10</sup><sup>−</sup>4; MDD <sup>=</sup> 10 km, <sup>α</sup> <sup>=</sup> <sup>2</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> HDD <sup>=</sup> 50 km, <sup>α</sup> <sup>=</sup> <sup>4</sup> <sup>×</sup> <sup>10</sup><sup>−</sup>5). The 'Species of Special Interest' function was used to increase the value of cells overlapping with *CAWA presence 2005–2009* and/or *CAWA presence 2010–2015.* Finally, we subtracted calculated *Model uncertainty* from all cells. After applying all features and functions, cells were prioritized for selection by the CAZ algorithm. LC2 included all features and functions from LC1. We then added *Protected areas* as a feature (weight 1.0) and used 'Matrix Connectivity' to prioritize selection of cells containing high Canada Warbler densities within 5 km of protected areas (distance specified by applying weighting factor = 0.1). The 'Matrix Connectivity' function multiplies the value of a cell based on its connectivity to other specified features, with features being considered connected if they fall within the specified dispersal distance [18].

LC3 included all features and functions from LC1 with the addition of the 'Ecological Interactions' function, which prioritizes areas based on connectivity between a pair of features [18]. We thus prioritized areas where high densities of Canada Warbler in both *CAWA climate baseline* (weight 0.75) and *CAWA climate 2050* (weight 0.75) were connected based on dispersal distance over 36 years (dispersal ability kernel is represented by <sup>β</sup> [18]; LDD <sup>=</sup> 5 km dispersal/year, 180 km total, <sup>β</sup> <sup>=</sup> 1.1 <sup>×</sup> <sup>10</sup><sup>−</sup>5; MDD <sup>=</sup> 10 km dispersal/year, 360 km total, <sup>β</sup> = 5.56 <sup>×</sup> 10−6; HDD = 50 km dispersal/year, 1800 km total, <sup>β</sup> = 1.1 <sup>×</sup> 10−6). We weighted climate features at 0.75 because assumptions about the future are less certain than present estimates (Moilanen and Arponen 2011). LC4 included all features and functions from LC2, with the addition of the 'Ecological Interactions' function used in LC3, in order to prioritize areas connecting high current and future populations of Canada Warbler with each other and with *Protected areas*.

FM1 included all features and functions from LC1 and added the 'Administrative Units' function, which recognizes that conservation decisions can be limited by administrative boundaries [18]. *Working lands* were given a weight of 0.5 to prioritize these areas for selection while maintaining connectivity with outside areas. The 'Matrix Connectivity' function was used to prioritize areas connecting high population density and *Wet-poor habitat.* FM2 included all features and functions from LC1 while omitting the 'Species of Special Interest' function to avoid prioritizing cells with Canada Warbler occurrences. 'Matrix Connectivity' was added to prioritize areas connecting high population density and *Upland habitat* (weight = 0.5). This function is based on the assumption that upland areas have desirable timber value and a greater potential to support Canada Warbler populations in 10–30 years after harvest if connected to dispersal sources and managed appropriately [13] (see Supplementary Materials). We added a second 'Matrix Connectivity' function to disincentivize prioritization of areas connecting high population density with *Wet-poor habitat* (weight = −0.5)*,* and a third 'Matrix Connectivity' function deterring prioritization of areas connecting high population density with *Protected areas* (weight = −0.5).

FM3 was not completed in Zonation, but rather was derived by subtracting the FM2 solution raster from the FM1 solution raster in ArcGIS. This was intended to identify areas for timber harvest with the lowest risk of harm to populations and maximum economic opportunity.

#### 2.3.3. Scenario Comparisons

For scenarios LC1–4 and FM1–2, we plotted performance curves representing conservation efficiency (the proportion of current predicted Canada Warbler population protected as a function of area protected, or not harvested in the case of FM2, at each dispersal distance estimate). Because of the different objectives of the scenarios, and functions applied, it was not sensible to compare conservation efficiency of all scenarios. We compared efficiency curves for scenarios with and without climate change but with all other factors being equal (LC2 and LC4). We also compared FM1 and FM2. Finally, we compared mean differences in cell-level rankings for all land conservation scenarios and estimated dispersal distances.

#### **3. Results**

For high resolution maps and data for all scenarios, visit https://github.com/borealbirds/cawa-bcr-14. Of the land conservation scenarios (Figure 2), areas that were consistently prioritized in both current and future climate scenarios included lands in central New Brunswick and the southeastern part of Québec along the border with Maine. All land conservation scenarios prioritized cells with recent and repeated observations of Canada Warbler, but these cells had a minimal effect on conservation efficiency or overall distribution of priority areas (e.g., Figure 3). Relatively few areas in Nova Scotia were prioritized in both current and future climate scenarios, except for small areas close to two large national parks. Adding connectivity to protected areas had little impact on the geographic distribution of areas prioritized in the current climate scenario (LC2), but had a more dramatic effect in the future climate scenario (LC4).

**Figure 2.** Zonation scenarios to prioritize areas for land conservation for Canada Warbler in Bird Conservation Region 14 at three natal dispersal distances: low dispersal distance (LDD, 5 km), medium dispersal distance (MDD, 10 km), and high dispersal distance (HDD, 50 km). Scenario descriptions given in Table 1. Priority rank for conservation scaled from blue (highest) to brown (lowest). Protected areas indicated by outlined polygons.

**Figure 3.** Zonation scenario to prioritize areas for land conservation with high current population density of Canada Warbler which are connected to protected areas at a medium natal dispersal distance (10 km). Crosshatched polygons indicate protected areas. Blue squares were highly prioritized because they included observations of Canada Warbler made during point count surveys between 2005 and 2015.

Increasing the dispersal distance estimate increased the aggregation of prioritized areas, with clusters of priority sites being larger and fewer in number for high dispersal distance scenarios, and smaller and more widely distributed for low dispersal distance scenarios (Figure 2). Overall, in the high dispersal distance scenarios, areas in Nova Scotia were rarely prioritized, with the largest aggregations of prioritized cells occurring in central New Brunswick, and to a lesser extent, southeastern Québec for current climate scenarios, and almost exclusively in Québec for future climate scenarios.

With the forest management scenarios (Figure 4), the emphasis was on forestry tenures; thus results are most appropriately applied by individual managers to compare the relative values of areas within those tenures when considering where to engage in responsible forest management practices (guidance is available in an accompanying technical report [33]) (see Supplementary Materials). Both FM1 and FM2 were affected by dispersal distance, shifting from many prioritized areas in Québec under LDD and MDD scenarios to prioritizing almost exclusively areas in New Brunswick under HDD scenarios.

When subtracting retention areas (FM1) from harvest priorities (FM2), the resulting scenario FM3 indicated that areas in the northern Gaspé Peninsula of Québec would be most effective at maximizing harvest value while minimizing potential loss of Canada Warbler. As with the land conservation scenarios, the prioritized areas were more aggregated with increased dispersal distance.

#### *Scenario Comparison*

For all scenarios, efficiency (proportion of current predicted Canada Warbler population retained per unit area) declined with higher dispersal distance estimates. Including climate change reduced land conservation scenario efficiency with respect to the current population (Figure 5). Efficiency was not directly comparable between forest management and land conservation scenarios due to differing objectives, as evidenced by the inverted efficiency curves for FM1 and FM2 (Figure 5), the former of which was designed to prioritize areas to avoid during forest management with the latter prioritizing areas to target for harvesting.

**Figure 4.** Zonation scenarios to prioritize areas for responsible forest management for Canada Warbler within forestry tenures in Bird Conservation Region 14 at three natal dispersal distances: 5 km (LDD), 10 km (MDD), and 50 km (HDD). Scenario descriptions given in Table 1. Priority rank for management scaled from blue (highest) to brown (lowest). Protected areas indicated by gray polygons.

**Figure 5.** Response curves comparing efficiency (proportion of Canada Warbler population density retained per unit area) for two forest management scenarios (left panel) and two land conservation scenarios (right panel). Scenario descriptions given in Table 1. Curves shown at three dispersal distances: 5 km (LDD, solid line), 10 km (MDD, hashed line), and 50 km (HDD, dotted line).

For land conservation scenarios, cell-level rankings, which range from 0 to 1, differed more between MDD and HDD scenarios (mean = 0.10, SD = 0.90) than between LDD and MDD scenarios (mean = 0.04, SD = 0.04), with greater variation in differences also found between MDD and HDD (Table 3).

**Table 3.** Percent change in mean cell-level rankings across Zonation prioritization scenarios to support land conservation for the Canada Warbler when comparing between different estimated dispersal distances: 5 km (LDD), 10 km (MDDs), and 50 km (HDD). Scenario descriptions given in Table 1.


#### **4. Discussion**

We generated land conservation and forest management prioritization maps for the Canada Warbler in the Canadian portion of the Atlantic Northern Forest. This exercise identified several consistently prioritized areas for a range of stewardship objectives. In particular, central New Brunswick and southern Québec emerged as important areas for conservation and management under both current and future climate scenarios and when considering a range of possible dispersal distances. In general, areas farther from the coast were more frequently prioritized for conservation. Including projected effects of climate change on potential population density had dramatic effects in the scenarios, leading to almost no areas prioritized for conservation in Nova Scotia, and shifting priority areas northward. Within forest management tenures, the northern Gaspé Peninsula of Québec was consistently identified as the area where forest harvesting activities may be most practical while avoiding impacts to current Canada Warbler populations. This is consistent with Sólymos et al. [44], who predicted a lower average Canada Warbler population density in the Gaspé Peninsula compared to the rest of the study area.

Given this species' high conservation concern, our single-species exercise has potential to aid the rapid implementation of conservation and recovery action. Our method was somewhat unique in that landscape prioritization exercises are typically used for the assessment of biodiversity at large scales considering many different species (e.g., [19,51]). For a given taxon, Zonation results that are produced using only survey data could be different from those using species distribution models [52] and including both as input features may offer benefits. Although in our case adding recent locations of Canada Warbler only had a small impact on the overall solution, important differences were apparent at the level of individual cells (1 km2). Including this element is critical for land managers making decisions on small scales, as they indicate land parcels with persistent occupancy by Canada Warblers. These areas can thus be targeted for permanent land conservation and avoided during forest harvesting.

Although our prioritization scenarios provide insight into possible locations to target conservation and management activities, the spatial resolution of the scenarios (1 km) is too coarse to identify specific habitat patches. While two finer-scale Canada Warbler species distribution models have been constructed in this region [10,53], their coverage does not include the entire study area. Prioritized areas should be regarded as suggestions that require more detailed investigation through comparison with satellite imagery and ground-truthing before being considered candidates for conservation or management activity. To maximize efficacy, our analysis should be repeated as finer-scale spatial data

become available to support local management planning and to account for the small territory sizes of the Canada Warbler (average 1 ha; [16,54]).

#### *4.1. Accounting for Uncertainty*

The predictions made by species distribution models, frequently used as input features in conservation prioritization exercises, are inherently more uncertain in under-sampled areas [55,56]. We accounted for this uncertainty by discounting densities by the standard deviation of mean population density, and thus were able to focus priorities on areas of lower scenario uncertainty [24,57].

We also attempted to account for uncertainty about dispersal distances within Canada Warbler populations by evaluating each scenario using three different dispersal distance estimates. These estimates were intended to influence results by dictating the extent to which priority conservation and management areas were required to be in close proximity to protected areas and known Canada Warbler locations, and to influence the allowed distance between current and future projected Canada Warbler distributions. However, we also used dispersal distance estimates within Zonation's 'Distribution Smoothing' function, which aggregates priority areas based on the assumption that fragmented solutions are undesirable [18], thereby yielding results that were increasingly spatially aggregated with larger dispersal distance estimates. This led to large differences in geographic priorities across scenarios, suggesting the need for careful consideration of habitat connectivity requirements for Canada Warbler and other species of conservation concern. Given that dispersal is a key parameter in spatial conservation prioritization [50], our findings highlight the importance of considering a range of dispersal assumptions. Future studies may benefit from multiple scenarios designed to isolate the different ways in which dispersal estimates can influence results.

The increased level of spatial aggregation in conservation priorities with higher dispersal estimates led to less efficient solutions in terms of the current predicted Canada Warbler population that would be conserved per unit area. However, conservation efficiency was influenced less by assumptions about dispersal distance than by assumptions about climate change effects on species' distributions, due to the large discrepancies between predicted current and future suitable habitats for boreal species [17]. Our results are consistent with those of Stralberg et al. [24], who found that incorporating potential avian responses to climate change reduced conservation efficiency for current songbird populations. This supports the idea that planning to incorporate climate change increases the size required for protected areas to adequately conserve species [17,58,59].

#### *4.2. Maintaining Viable Populations on the Landscape*

Connectivity was important in identifying management opportunities through this prioritization effort because Canada Warblers cluster in multi-territory "neighborhoods" [16,60]. Canada Warblers use forest stands post-harvest more often if they are within 100 m of unharvested stands with conspecific breeders; in one study, the presence of conspecifics was found to be more important than habitat condition for predicting stand use [16]. Therefore, it is important to conserve or manage for areas of sufficient size to support a neighborhood, although the ideal size and configuration of such habitat patches is not yet known. It is also not known whether dispersal within and between such neighborhoods is important for Canada Warbler population viability, nor what barriers to functional connectivity [61] exist for this species in this region. However, observed population declines combined with stable breeding productivity [62] suggest that eastern populations may not be limited by quality or connectivity of breeding habitat.

Our scenarios that assumed a high dispersal distance (HDD) prioritized large areas in the northwest portion of the study region. In contrast, the low dispersal distance (LDD) scenarios prioritized more locations of smaller size across the Atlantic Northern Forest. Relying on the results from HDD scenarios alone could undermine the goals of individual provinces to maintain native species by favoring the conservation of fewer large populations rather than a larger number of small, spatially distributed populations. The latter may be preferable from the standpoint of maintaining genetic variability across the landscape [63]. Furthermore, our scenarios that considered future climate projections (particularly the HDD variants) assigned low priority to southern populations in Nova Scotia and southern Québec—areas that become less hospitable to Canada Warbler in a warmer climate. Thus, the LDD and current climate scenarios represent more conservative assumptions for conservation and management. This may be more appropriate for a species that is experiencing population declines, especially given range-wide projected increases in habitat suitability under climate change [17].

To guarantee long-term persistence of high-quality breeding areas for the Canada Warbler, information on the capacity of habitats to support viable populations through detailed spatial population viability analyses is critical [64]. Although Zonation uses connectivity of populations as a surrogate for viability [20], this may not be as relevant for passerines, who have greater mobility than taxa with small dispersal distances such as small mammals, plants, or colonial birds. By incorporating measures of connectivity and including recent and repeated observations of Canada Warblers to locate high priority areas, we may have been able to better prioritize the landscape for conservation of high-value habitat than what was accomplished by using species distribution models or survey data alone. We were not able to include data on habitat selection and reproductive success that may have given a more direct indicator of population viability, which is being used in ongoing studies (Burns & Reitsma, in revision; Amelie Roberto-Charron, pers. comm.; Junior Tremblay, pers. comm.). Future prioritization efforts should include results of assessments of reproductive success and population viability wherever possible in order to most accurately target areas to maintain population density on the landscape. Such data would be particularly valuable if comparing harvested and unharvested areas, to test the hypothesis of whether areas undergoing forest management represent ecological traps for this species [65].

The present study advances the understanding of conservation issues and management opportunities for Canada Warbler at a regional scale. Our results help to define priority areas for Canada Warbler land conservation and forest management, and in conjunction with the habitat guidelines for the Canadian portion of the Atlantic Northern Forest [13] (see Supplementary Materials), they provide a toolkit for managers to immediately locate areas for implementing conservation and management actions. We suggest that this approach, designed to support management objectives for a single species, be applied to other species. We particularly encourage managers to apply this prioritization approach to Canada Warbler populations in other BCRs in Canada and the U.S. to support persistence of the entire breeding population.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/1424-2818/12/2/61/s1, file 1: Westwood et al—2017—Guidelines for managing Canada Warbler habitat.pdf; file 2: Westwood, Reitsma, Lambert—2017—Prioritizing areas for Canada Warbler conservation and management.pdf.

**Author Contributions:** Conceptualization, J.D.L., A.R.W., and L.R.R.; methodology—development, spatial analysis, and statistics, A.R.W. and D.S.; writing—original draft preparation, review and editing, A.R.W., J.D.L., L.R.R., and D.S.; funding acquisition, J.D.L. and A.R.W. All authors have read and agree to the published version of the manuscript.

**Funding:** This research was funded by Environment and Climate Change Canada (ECCC), the U.S. Fish and Wildlife Service (USFWS), and High Branch Conservation Services Ltd. Article Processing Charges were funded by [TBA].

**Acknowledgments:** We acknowledge our funders for their support through various components of this project. Thank you to Alaine Camfield and Sean LeMoine of ECCC who were highly involved in the preparation of associated technical reports. The associated technical reports acknowledge the names of the many professionals across sectors (government, industry, academia, and Indigenous communities) who responded to consultations and surveys associated with this project, and for whose time and expertise we are grateful. We acknowledge the BAM Project's members, avian and biophysical data partners, and funding agencies, listed in full at www.borealbirds.ca/index.php/acknowledgements. The BAM project's database includes Breeding Bird Survey data and provincial Breeding Bird Atlas data, and we acknowledge the hundreds of skilled volunteers in Canada who have participated in these surveys over the years, and those who have served as provincial or territorial coordinators for the BBS. We thank the subject editor and two reviewers for a previous submission to the journal Avian Conservation and Ecology for their constructive comments, which refined the manuscript and led us to find a more appropriate home for applied conservation work. We are grateful to BAM members Péter Sólymos, and Trish Fontaine for provision of data and Nicole Barker and Samuel Haché for guidance and feedback. Finally, we thank Kara Pearson for GIS database assembly and Charlotte Harding for contributions to the U.S. habitat management guidelines.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Occupancy of the American Three-Toed Woodpecker in a Heavily-Managed Boreal Forest of Eastern Canada**

**Vincent Lamarre <sup>1</sup> and Junior A. Tremblay 1,2,\***


**Abstract:** The southern extent of the boreal forest in North America has experienced intensive human disturbance in recent decades. Among these, forest harvesting leads to the substantial loss of late-successional stands that include key habitat attributes for several avian species. The American Three-toed Woodpecker, *Picoides dorsalis*, is associated with continuous old spruce forests in the eastern part of its range. In this study, we assessed the influence of habitat characteristics at different scales on the occupancy of American Three-toed Woodpecker in a heavily-managed boreal landscape of northeastern Canada, and we inferred species occupancy at the regional scale. We conducted 185 playback stations over two breeding seasons and modelled the occupancy of the species while taking into account the probability of detection. American Three-toed Woodpecker occupancy was lower in stands with large areas recently clear-cut, and higher in landscapes with large extents of old-growth forest dominated by black spruce. At the regional scale, areas with high probability of occupancy were scarce and mostly within protected areas. Habitat requirements of the American Three-toed Woodpecker during the breeding season, coupled with overall low occupancy rate in our study area, challenge its long-term sustainability in such heavily managed landscapes. Additionally, the scarcity of areas of high probability of occupancy in the region suggests that the ecological role of old forest outside protected areas could be compromised.

**Keywords:** boreal forest; clear-cutting; conservation; forest management; old-growth forest; *Picoides dorsalis*; protected areas

#### **1. Introduction**

The boreal forest represents about 48% of the world's forested biomes [1]. Boreal landscapes are shaped by natural disturbances such as windthrows, forest fires and insect outbreaks that generate a complex mosaic in vegetation structure and composition [2,3]. For several decades, however, the exploitation of natural resources has been adding to natural disturbances in this biome [4] and has modified the structure and composition of boreal ecosystems, including the reduction of late-successional stands called old-growth forests [5,6].

Old-growth boreal forests are characterized, among several attributes, by irregular vertical and horizontal structures, large volumes of deadwood either standing (snags) or fallen (coarse woody debris) and in different decaying stages [7]. These attributes are considered key for hundreds of species that depend upon dead or decaying woody material during some part of their life cycle and that are found disproportionately in old-growth forests [8,9]. The temporal continuity of these ecological attributes is also an essential characteristic for many species [10,11]. For instance, a recent study highlights the continuous supply of large slightly decayed snags in specific old-growth forest type as a key element to provide temporal stability in the foraging habitat of the Black-backed Woodpecker, *Picoides arcticus* [12].

**Citation:** Lamarre, V.; Tremblay, J.A. Occupancy of the American Three-Toed Woodpecker in a Heavily-Managed Boreal Forest of Eastern Canada. *Diversity* **2021**, *13*, 35. https://doi.org/10.3390/ d13010035

Received: 17 December 2020 Accepted: 15 January 2021 Published: 19 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Deadwood provides foraging or breeding substrate for several vertebrate species [13]. Boreal woodpeckers, for example, are "ecosystem engineers" that create nesting cavities for other vertebrates and are considered indicator species for deadwood-associated biodiversity [14,15]. Among woodpecker species found within the boreal biome in North America, the American Three-toed Woodpecker, *Picoides dorsalis*, is the most strongly associated with continuous old spruce forests at the landscape scale [16–18]. Harvesting, especially of oldgrowth coniferous forests, is thought to be detrimental to the species by limiting the supply of deadwood [16,18–20]. However, most of the studies on the species in eastern Canada occurred in regions where forests were harvested for the first time and where mature and old stands were still relatively abundant amongst residual forests. In heavily-managed forest landscapes, the persistence of the species remains uncertain. Indeed, in a recent attempt to gain insight on the breeding ecology of American Three-toed Woodpecker in heavily managed forest at the southern edge of its breeding range (New-Brunswick, Canada), Craig et al. [21] reported the species in only 5.9% of the playback stations, and found no predictors of site occupancy although most nests were found in recently dead black spruce trees.

Here, we assess the influence of habitat characteristics at the stand and landscape scales on the occupancy of American Three-toed Woodpecker during the breeding period in a heavily-managed, unburned boreal forest landscape of eastern Canada, and infer the probability of occupancy of the species at the regional scale. We expect that American Three-toed Woodpecker occupancy would be favoured by old spruce forest stands but negatively affected by recently harvested forest stands and higher forest fragmentation, and that the probability of occupancy of the species at the regional scale would be low.

#### **2. Materials and Methods**

#### *2.1. Study Area*

This study was conducted at the Forêt Montmorency (Université Laval's Experimental Forest), and the Réserve Faunique des Laurentides, Québec, Canada (47◦4 N. 71◦0 W.; Figure 1). The study area covers approximately 210 km<sup>2</sup> within the balsam fir-white birch bioclimatic domain in the continuous boreal forest subzone [22], which represents the southernmost section of the boreal forest in eastern Canada [23]. More precisely, the study area is classified within the high-elevation balsam fir–white birch zone with elevation ranging between 600 and 1100 m [24]. Forest stands within the study area are dominated by balsam fir *Abies balsamea* and black spruce *Picea mariana*, the latter being more abundant at higher elevation. Companion species are white birch (*Betula papyrifera)*, white spruce (*Picea glauca),* and tamarack *(Larix laricina)*. In the region, average age-class structure of the natural variability observed over the last few centuries includes 86% of old forest stands [25]. However, intensive logging during the 20th and early 21st century led to a net reduction in the cover of mature and old forest stands [26]. During the 1909–2005 time frame, 105% of the high-elevation boreal landscape has been harvested with some areas logged twice [27]. Consequently, the studied area is now dominated by young regenerating stands (0–19 years; 39%), while the cover of old forests (>90 years) is approximately 25%. Young closed (20–59 years) and mature (60–89 years) forests as well as non-forested lands (e.g., lakes, wetlands, etc.) respectively cover 11, 16 and 9% of the study area. Dominant cover tree species of old forest stands are balsam fir (55%), black spruce (42%) and tamarack (3%).

**Figure 1.** Habitat composition and location of surveyed stations where American Three-toedWoodpecker was detected (black circles) and not detected (grey circles) during occupancy surveys conducted between 2016 and 2017 within the Forêt Montmorency and the Réserve Faunique des Laurentides (Québec, Canada). Boundaries of the boreal zone [28] are delimited by the green zone in the left insert. The black line delineates the boundaries of the study area in the large insert.

#### *2.2. Woodpecker Surveys*

In 2016 and 2017, we conducted woodpecker occupancy surveys during American Three-toed Woodpecker breeding season at stationary stations distributed along forest roads. We did not stratify the sampling *per se* but rather targeted areas with clusters of unharvested patches within the vicinity of old forest stands. Surveyed stations were distanced by a minimum of 255 m (mean ± standard deviation: 554 ± 442 m; range 255–3414 m), while mean elevation at stations was 889 ± 37 m (range 819–1015 m). Surveys lasted five minutes (300 s), during which conspecific playbacks were displayed, and we noted the time elapsed since the start of the playback when an individual was observed or heard. Similar to the method used by Craig et al. [21], the 300-s conspecific display was divided into two 120-s observation periods which constituted the two visits used in occupancy analyses (see below). Both observation periods were interspersed by a 60-s pause. Playbacks were mixed from recordings obtained from xeno-canto [29] and the Macaulay Library [30]. Playbacks could be heard to approximately 450 m by the human ear in the field. We conducted surveys between 5:47 a.m. and 3:32 p.m., when precipitations were absent or minimal and wind speed ≤ 3 on the Beaufort scale (12–19 km h<sup>−</sup>1). Afterward, we accessed three-hour mean wind speed from the nearest weather station located 18 km northwest of the study area for further analysis [31].

#### *2.3. Habitat Characteristics*

We investigated the influence of habitat characteristics that were likely to affect the occupancy of American Three-toed Woodpecker by calculating vegetation covariates within buffers of 250 and 750 m-radii centered on each surveyed station, representing respectively the stand and landscape scales. The landscape scale was selected based on home range size estimates obtained for one female and two male American Three-toed Woodpeckers during the nesting period in the study area, which averaged 177 ha (J.A. Tremblay, unpublished data). Vegetation covariates were calculated using the eco-forestry layers available as of 2017 for the study area [32]. At the stand scale, we calculated the area (ha) covered by recent clear-cut stands *clear\_cut*. At the stand and landscape scales, we measured the area (ha) covered by old forest (≥90 years) dominated by black spruce *old\_spruce* and mean forest age *mean\_age*. We weighted *mean\_age* by the area covered by even-aged stands within the buffers around each surveyed station. We calculated the standard deviation of mean stand age *sd\_age* as an index of habitat fragmentation at the landscape scale. Based on the literature in eastern Canada [16–18], the area covered by non-forested lands in the landscape, the area covered by young (20–59 years) and mature (60–89 years) forests at the stand and landscape scales have been considered marginal for the American Threetoed Woodpecker in the selection of its habitat during the breeding season and thus were excluded from occupancy analyses (see below). We nevertheless compared mean value of all habitat characteristics between stations where the American Three-toed Woodpecker was detected, not detected and a set of 185 random stations distributed randomly into the study (Table 1) area using a Kruskal–Wallis test in R statistical environment version 4.0.2 [33].

**Table 1.** Description and mean values ± standard deviation (s.d.) of habitat characteristics measured at the stand (250 m-radii) and landscape (750 m-radii) scales at stations where the American Three-toed Woodpecker was detected, not detected and at random stations during occupancy surveys in managed boreal forest in eastern Canada. Shared letters indicate no significant differences in habitat characteristics among stations (Kruskal–Wallis test). Habitat characteristics in bold are included in occupancy analyses.


#### *2.4. Occupancy Analyses*

We used single-species occupancy modeling in the *unmarked* R library to investigate the influence of habitat characteristics on the probability of occupancy of American Three-toed Woodpecker [34,35]. We converted survey detections of American Three-toed Woodpecker into presence–absence data. To account for imperfect detection probability during surveys, we first estimated the effect of year, Julian day, wind speed (m s−1) and hour of the day on detection probability of American Three-toed Woodpecker. We used a two-step approach [36] to determine which detection parameter(s) to retain in the occupancy models. We first estimated the effect of each detection parameter by ranking univariate models in which occupancy was held constant. We also included a null model in which detection and occupancy were held constant (Table A1). The five candidate models were ranked based

on the second-order Akaike's information criterion AICc [37]; using the aictab function of the *AICcmodavg* R library [38].

We developed a set of 10 biologically relevant candidate occupancy models representing hypotheses to investigate the influence of habitat characteristics on the probability of occupancy of American Three-toed Woodpecker (Table 2). This set also included a null model of constant occupancy. Strongly correlated covariates (|*r*| ≥ 0.60; Table A2) were not included in the same candidate model to reduce multicollinearity. We ranked candidate models based on the AICc and we made inference on the top models with ΔAICc < 2. Parameter-averaged predictions of covariates appearing in the most parsimonious models were calculated over their measured range while holding other variables at their mean value. Finally, we inferred a predicted probability of occupancy of American Three-toed Woodpecker from the top models (Table 2) at the regional scale (30 km × 30 km square zone) within the balsam fir–white birch bioclimatic domain.

**Table 2.** Candidates models, number of parameters (*k*), second-order Akaike's information criterion (AICc), ΔAICc, Akaike weights (*ω*), log-likelihood (LL) of the candidate models assessing occupancy (ψ) and detection probability (*p*) of American Three-toed Woodpecker in a managed boreal forest in eastern Canada.


#### **3. Results**

Between 19 May and 13 July 2016 and between 26 May and 29 June 2017, we broadcasted American Three-toed Woodpecker playbacks at 185 stations, including 35 revisited stations in 2017. American Three-toed Woodpecker was detected at 22 (11.9%) of the 185 stations. Forest composition and structure differed between sampled and random stations since we targeted our sampling effort towards clusters of older residual forest patches distributed close to non-forested lands such as waterbodies and wetlands in the study area (Table 1). Indeed, mean forest age at the stand and landscape scale and the area covered by old spruce at the stand scale were lower at random than sampled stations. In addition, the area covered by mature stands at both scales was also lower at random stations, while the dominance of young stands tended to be greater. Finally, non-forested habitats in the landscape were less important at random than sampled stations (Table 1).

Time of the day was retained in the most parsimonious model and accounted for 42% of model selection weight (Table A1). The models including year (ΔAICc = 1.61) and date (ΔAICc = 1.95) were equivalent and did not differ from the null model (ΔAICc = 1.77). Only time of the day influenced the probability of detection of American Three-toed Woodpecker (*p*time: 0.37, 95% C.I.: [0.03, 0.71]) and therefore was the only detection parameter included in occupancy models. Mean occupancy rate of American Three-toed Woodpecker in the study site was 0.17 ± 0.04, and mean detection probability when the time of the day was fixed to its mean value was 0.48 ± 0.11.

Two models had a substantial level of empirical support in influencing occupancy of the American Three-toed Woodpecker (ΔAICc < 2; Table 2). Together, these models accounted for 68% of the AIC weight. At the stand scale, occupancy of American Three-toed Woodpecker decreased with increasing recent clear-cut area (ψclear\_cut\_250: −0.20, 95% C.I.: [−0.34, −0.06], Figure 2a). The area covered by recent clear-cuts at the stand scale was also lower at stations where American Three-toed Woodpecker was detected compared to stations where the species was not detected or random stations (Table 1). At the landscape scale, occupancy of the species tended to increase with an increase in the area covered by old forest dominated by black spruce (ψold\_spruce\_750: 0.04, 95% C.I.: [0.00–0.08], Figure 2b), with increasing uncertainty in confidence intervals > 40 ha most likely due to the scarcity of large old spruce forest stands in our study area.

**Figure 2.** Influence on the occupancy of American Three-toed Woodpecker in managed boreal forest in eastern Canada of the area (ha) covered by (**a**) recent clear-cut stands at the stand scale (250 m radii) and (**b**) old spruce stands at the landscape scale (750 m radii). The shaded grey area represents 95% confidence interval, and symbols represent the distribution of raw data with symbol size proportional to the (log + 1) number of observations.

Extrapolating results inferred from our top occupancy models (Table 2) to a larger extent within the balsam fir–white birch bioclimatic domain showed that only 12.3% of the region had moderate- to high-predicted probability of occupancy (>0.5) of the American Three-toed Woodpecker. About half of these areas (6.7%) were located within protected areas. Areas with high-predicted probability of occupancy (>0.75) of the species represented 4.7% of the region where only 1.0% were in managed forests outside protected areas (Figure 3).

**Figure 3.** Predicted probability of occupancy of American Three-toed Woodpecker at the regional scale within the balsam fir–white birch bioclimatic domain. Probability of occupancy is inferred from the top-ranking models (ΔAICc < 2; Table 2). The dashed and solid black lines respectively delineate the boundaries of the study area and protected areas and white polygons indicate non-forested lands.

#### **4. Discussion**

In a heavily-managed landscape at the southern edge of the boreal forest, the American Three-toed Woodpecker exhibited a low mean occupancy rate (0.17 ± 0.04) where its probability of occupancy decreased rapidly with increasing area of recent clear-cut at the stand scale. The amount of old spruce forest was positively associated with species occupancy at the landscape scale. Our results suggest that logging history in the region may have created an unsuitable forest landscape (i.e., dominated by younger age-classes and fragmented residual older forest stands) to support a long-term population of American Three-toed Woodpecker.

Habitat associations of the American Three-toed Woodpecker vary across its range. In a meta-analysis focusing on successional trajectories of bird communities following fire and harvesting in boreal forests of western Canada, Schiek and Song [39] report the American Three-toed Woodpecker more common in older mixed wood and white spruce. In western Quebec, Imbeau and Desrochers [16] document a positive association between the American Three-toed Woodpecker with the amount of old spruce forest in continuous coniferous forest. Similarly, Cadieux and Drapeau [18] report a higher occurrence of the species in black spruce forest stands older than 90 years and no occurrence in younger coniferous stands. Accordingly, we find that the occupancy of the American Three-toed Woodpecker during the breeding period increases with old spruce forest stands in a managed landscape. This result agrees with our prediction as it has been shown that the

American Three-toed Woodpecker preferentially forages on relatively large and senescent or recently dead coniferous trees to access bark-associated beetles (Scolytinae) by scaling layers of bark with strong preference for black spruce in eastern Canada [17,40,41].

Recent clear-cut at the stand scale have a negative effect on the occupancy of the American Three-toed Woodpecker. This result agrees with our expectation and with previous studies about the sensitivity of the species to forest harvesting, especially in the eastern part of its distribution range [19,20,41,42]. Accordingly, in the black spruce moss domain of eastern Canada, recent clear-cuts have a lower density of snags, and the American Three-toed Woodpecker is solely found in old-growth forest [19]. In addition to a reduction in the deadwood abundance, long-term recruitment of large-diameter snags is compromised in clear-cut stands [43]. In a closely-related species, the Black-backed Woodpecker, old coniferous forests with a greater volume of recently dead trees than adjacent recent cuts are selected for foraging during the breeding period [44]. Similar foraging avoidance of recent clear-cut areas may be occurring in our study site for the American Three-toed Woodpecker. Indeed, timber harvesting peaked during the period 1996–2004 and ended in 2009, and although we did not quantify the volume and decay stages of deadwood, it is reasonable to think that most deadwood in clear-cuts had entered late decay classes and was of low quality for foraging American Three-toedWoodpecker at the time of our study.

Our results do not report a predictive effect of habitat fragmentation at the landscape scale on the occupancy of the American Three-toed Woodpecker. Forest edge avoidance by foraging American Three-toed Woodpecker has been reported, where high-quality substrates near stand edges are used less frequently than available [41]. In addition, movements of foraging woodpeckers also appear to be constrained in residual forests following harvesting [16,41]. Hence, the effect of habitat fragmentation seems to act on a finer scale, and our results suggest that at a larger scale (i.e., landscape), the amount of old forests may strongly influence the occurrence of the species. For instance, most of our detection of the species occurred close to forest patches of mature or old forests rather than residual strips of forest (Figure 1).

Only 1% of the forest stands in managed forests at the regional scale show high predicted probability of occupancy of the American Three-toed Woodpecker. This is likely a consequence of forestry practices, mainly driven by clear-cutting during the last century which has led to a substantial reduction of old forest cover in the region [26,27]. Such practices can hardly sustain biodiversity associated with old coniferous forest stands, especially for species with large home ranges. For example, habitat alteration from anthropogenic activities is a threat of high concern for populations of Woodland Caribou *Rangifer tarandus caribou* across Canada's boreal biomes, including the local population in our study area, which is considered "not self-sustaining" [45]. Fennoscandia previously experienced a similar situation where old-growth forests have almost completely disappeared [46,47], and it is estimated that 30–50% of the red-listed species in these regions are associated with old-growth forest attributes [48,49]. In boreal ecosystems, conservation strategies cannot be based solely on a network of protected areas but rather on how unprotected areas are managed [50]. Hence, management practices mimicking the attributes of old-growth forests by ensuring a continuous recruitment of deadwood appear important for maintaining species associated with old spruce forests. Within these practices, partial harvesting may be efficient in maintaining a relatively high abundance of deadwood and associated deadwood-dependent species such as the American Three-toed Woodpecker [51,52], although the efficiency of such practices has not yet been assessed. The ecological role of old forests in our study area seems to be altered, and the situation may require revised management practices and a passive restoration of the ecological integrity of old forests outside protected areas in the region *sensu* [53].

#### **5. Conclusions**

With 185 stations sampled over two breeding seasons, our study on the occupancy of the American Three-toed Woodpecker is one of the first conducted in a heavily harvested landscape at the southern edge of the boreal forest. The overall low occupancy of the species within our study area raises questions about the long-term sustainability of such heavily managed landscapes for the American Three-toed Woodpecker. Areas of highpredicted occupancy of the species in the studied region are mostly found in protected areas, providing evidence that heavy forest harvesting is a detrimental driver that likely contributed to the significant long-term declining trends of the species over the 1970–2019 period in the province of Québec (−0.96%.year−1; 95% C.I.: [−0.70–−1.22]; [54]) and at a larger scale in the boreal hardwood transition bird conservation region (−3.5%.year−1; 95% C.I.: [−2.0–−5.2]; [55]). The southern part of the boreal forest in eastern Canada, where we conducted our study, has experienced one of the most important human disturbances in the past decades, mainly related to forest harvesting [4]. We pledge for more detailed studies on this discrete keystone species in different regions with varying forest harvesting intensity, at the stand scale (going from clear-cutting to partial harvesting) and landscape scale (from pristine to heavily harvested), focusing on demographic parameters (i.e., reproductive success and survival) and habitat selection.

**Author Contributions:** Conceptualization, J.A.T.; methodology, J.A.T.; formal analysis, J.A.T., V.L.; investigation, J.A.T., V.L.; resources, J.A.T.; data curation, V.L.; writing—original draft preparation, V.L.; writing—review and editing, J.A.T., V.L.; visualization, V.L.; supervision, J.A.T.; project administration, J.A.T.; funding acquisition, J.A.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Environment and Climate Change Canada.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are openly available in https://doi. org/10.5061/dryad.9kd51c5gc.

**Acknowledgments:** We are grateful to André Desrochers for sharing the code to calculate habitat characteristics and to Mark J. Mazerolle for advice related to statistical analyses. We are also thankful to the personnel at the Forêt Montmorency for their support and to Michel Robert for the loan of allterrain vehicles. Finally, we thank Francis Lessard and Julien St-Amand for their help with fieldwork.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A**

**Table A1.** Model selection based on the AICc for the estimation of detection probability (*p*) of American Three-toed Woodpecker in managed boreal forest in eastern Canada. The number of parameters (*k*), second-order Akaike's information criterion (AICc), ΔAICc, Akaike weights (*ω*), log-likelihood (LL) are included in the table and occupancy (ψ) is held constant for each candidate model.



**Table A2.** Correlation among habitat characteristics included in occupancy analyses.

#### **References**


## *Article* **Implications of Historical and Contemporary Processes on Genetic Differentiation of a Declining Boreal Songbird: The Rusty Blackbird**

**Robert E. Wilson 1,\*, Steven M. Matsuoka 1, Luke L. Powell 2,3,4,5, James A. Johnson 6, Dean W. Demarest 7, Diana Stralberg 8,9 and Sarah A. Sonsthagen <sup>1</sup>**


**Abstract:** The arrangement of habitat features via historical or contemporary events can strongly influence genomic and demographic connectivity, and in turn affect levels of genetic diversity and resilience of populations to environmental perturbation. The rusty blackbird (*Euphagus carolinus*) is a forested wetland habitat specialist whose population size has declined sharply (78%) over recent decades. The species breeds across the expansive North American boreal forest region, which contains a mosaic of habitat conditions resulting from active natural disturbance regimes and glacial history. We used landscape genomics to evaluate how past and present landscape features have shaped patterns of genetic diversity and connectivity across the species' breeding range. Based on reducedrepresentation genomic and mitochondrial DNA, genetic structure followed four broad patterns influenced by both historical and contemporary forces: (1) an east–west partition consistent with vicariance during the last glacial maximum; (2) a potential secondary contact zone between eastern and western lineages at James Bay, Ontario; (3) insular differentiation of birds on Newfoundland; and (4) restricted regional gene flow among locales within western and eastern North America. The presence of genomic structure and therefore restricted dispersal among populations may limit the species' capacity to respond to rapid environmental change.

**Keywords:** *Euphagus carolinus*; genetic diversity; boreal; glacial refugia; phylogeography

#### **1. Introduction**

The spatial organization of suitable habitat across the landscape via historical and contemporary events plays an important role in the maintenance of both genetic and demographic connectivity within plant or animal populations. Discontinuities in habitat, whether from physical barriers or from natural or anthropogenic disturbances, fragment populations into smaller units and can diminish levels of connectivity when (1) the distance

**Citation:** Wilson, R.E.; Matsuoka, S.M.; Powell, L.L.; Johnson, J.A.; Demarest, D.W.; Stralberg, D.; Sonsthagen, S.A. Implications of Historical and Contemporary Processes on Genetic Differentiation of a Declining Boreal Songbird: The Rusty Blackbird. *Diversity* **2021**, *13*, 103. https://doi.org/10.3390/ d13030103

Academic Editor: Stacy A. McNulty

Received: 9 February 2021 Accepted: 16 February 2021 Published: 25 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

between suitable habitat patches exceeds the dispersal capacity of individuals or (2) dispersing individuals do not successfully reproduce. In this way, spatial habitat heterogeneity influences effective dispersal (dispersal followed by reproduction) by individuals, which in turn impacts gene flow and population dynamics (i.e., potential outcome of dispersal, [1,2]). Isolation, whether by distance or environment, often results in lower levels of intra-population genetic variation and higher levels of inter-population genetic differentiation across a species' range, which can have major short- and long-term implications on population persistence. For example, isolated populations with reduced genetic diversity may lose adaptive potential, accumulate deleterious mutations, or experience increased inbreeding [3–7]—all of which can reduce individual fitness and increase vulnerability of populations to decline and extirpation [8–10].

While habitat configuration can have a profound influence on population structure and dynamics, species-specific responses to habitat heterogeneity vary widely due to differences in life history characteristics [11]. In some species, reduced connectivity among habitats may be offset by either high dispersal capabilities [12,13] or the presence of relatively continuous habitat during critical parts of the annual cycle that promote genetic exchange (e.g., mating). Even in highly vagile migratory species, however, habitat fragmentation and availability can contribute to genetic differentiation [14,15]. Thus, the capacity to move may be insufficient alone to offset the negative effects of habitat loss. Effective dispersal (i.e., with reproduction) can have a stronger effect on population demographics than the dispersal event alone [16–18]. Because of this, conservation strategies may need to promote site conditions favorable for successful reproduction in addition to general landscape connectivity, as gene flow is essential for population persistence in a changing environment [16–18]. However, determining the effective outcome of dispersal events is often difficult or nearly impossible with banding or telemetry data alone as neither can directly infer successful reproduction. Genetic signatures can provide insights into the success of natal dispersal and connectivity among breeding areas, which are relevant to the conservation of populations. Population genomics can help identify areas where connectivity across the landscape enriches genetic diversity and enhances the resiliency of populations to environmental perturbation. It can also identify where contemporary or historical limitations in dispersal have led to genomic structuring and distinct populations that may require specific management strategies to remain viable. Thus, population genomics provides a powerful approach to understanding implications of dispersal and can fill information gaps in traditional movement data, especially in migratory birds that nest in remote regions (see references within [15]) such as the vast boreal forest biome of North America.

North America's boreal forest biome contains ≥ 25% of the world's wetlands and intact forests [19,20], which provide important habitats to over 300 bird species during the breeding season [21]. While migratory songbirds across North American have undergone concerning population declines, boreal nesting species have exhibited some of the most dramatic declines over the past half century [22,23]. Rusty blackbirds (*Euphagus carolinus*) are an unfortunate example of this pattern; since 1966 the global population size is estimated to have decreased by 78% [23,24], with the decline likely ongoing since the late 19th century [25,26]. Loss of wooded wetlands on the wintering grounds in the southern United States is suspected as a principal driver of these declines [27]. However, methylmercury contamination on the breeding grounds [28–30], conversion of wetland habitats used during migration [31], alteration of boreal wetland breeding habitats, and climate change are also likely contributing factors [27,32,33].

Little is known regarding patterns of genomic connectivity and differentiation among nesting areas of rusty blackbirds, owing to their expansive and remote distribution spanning the boreal forest biome from Alaska eastward to Newfoundland and northern New England. However, banding, migration tracking, and stable-isotope analyses indicate a general migratory divide. Rusty blackbirds nesting in Alaska and western Canada migrate west of the Appalachian Mountains to winter in the Mississippi Alluvial Valley, while

birds nesting in eastern Canada and the northeastern United States migrate east of the Appalachian Mountains and winter on the Southeastern Coastal Plain [34–36]. There have been no phylogenetic studies on the species to date, but birds breeding on Newfoundland and Magdelen Islands, Quebec and wintering in South Carolina have been described as a putative subspecies (*E*. *c*. *nigrans*) that is phenotypically distinct from the nominate form breeding over the remainder of the range [37].

Here we present a landscape genomic approach to examine the influence of historical (Last Glacial Maximum, LGM) and contemporary processes on patterns of diversity within the rusty blackbird, using reduced representation genomic (double-digest restriction-site associated DNA sequences, ddRAD; biparentally inherited) and mitochondrial DNA data (mtDNA, maternally inherited). Much of the present-day geographic extent of the Nearctic boreal biome was covered by ice sheets during the LGM ~20,000 years ago. Consequently, the post-glacial colonization of the region by flora and fauna expanding their ranges from south of ice sheets or from ice-free glacial refugia such as Alaska (Beringia) and the Atlantic Shelf [38] has influenced how genetic diversity within species is arrayed across boreal landscapes (e.g., [39–42]). Within glacial refugia, populations that diverged in isolation during the LGM generally harbor greater levels of genetic diversity than populations in areas recolonized after glacial retreat, except in areas where lineages from separate refugia intermixed [43]. As the present-day distribution of rusty blackbirds spans the entire North American boreal biome including previously glaciated and unglaciated areas, we hypothesize that (1) rusty blackbirds likely contracted to at least two refugia during the LGM as suggested for other avifauna (e.g., [44]). Further, we hypothesize that (2) the putative subspecies of rusty blackbird (*E*. *c*. *nigrans*, [37]) residing on the island of Newfoundland would be differentiated from eastern populations on the mainland through insular isolation, LGM isolation in the Atlantic Shelf refugium, or isolation through other physical barriers (e.g., [45,46]).

Just as the repeated glacial and interglacial cycles of the 2.5 million-year Pleistocene epoch were a powerful force shaping patterns of genetic diversity in boreal forest birds [47], contemporary genetic diversity and structure are also dependent on behavioral and biological aspects of individual species [48] that influence their capacity to move across the landscape [49,50]. Indeed, interpopulation variation in dispersal traits is often related to landscape structure (see [51]). The boreal biome is a mosaic of wetland complexes, upland forests, and montane areas and therefore comprises a naturally fragmented landscape for wetland habitat specialists, such as the rusty blackbird [52]. Although rusty blackbirds are able to move long distances and traverse large gaps in suitable habitat during migration (e.g., Lake Erie, [31], western cordillera, [53]), the broad level migratory connectivity observed [35] may equate to some degree of philopatry within each region (i.e., return to natal area to reproduce). We therefore hypothesize that (3) rusty blackbirds will display genetic structure within western and eastern North America.

#### **2. Materials and Methods**

#### *2.1. Sampling and DNA Extraction*

Blood was sampled from 205 rusty blackbirds captured from 5 May to 15 July 2009– 2018 in mist nets placed near their nests or in post-breeding foraging areas across their breeding range (Figure 1, see Table 1 for sample sizes and locality names). Genomic DNA was extracted using a DNeasy Blood and Tissue kit following the manufacturer's protocols (Qiagen, Valencia, CA, USA). Extractions were quantified using a Broad Range Quant-iT dsDNA Assay Kit (Thermo Fisher Scientific, Inc., Waltham, MA USA).

**Figure 1.** Nesting distribution of rusty blackbirds with sampling locations color-coded and numbered (see Table 1). We note that location 13 (white circle) includes sample locations in Vermont and New Hampshire, USA, and location 12 (grey circle) includes both Nova Scotia and New Brunswick, Canada. The scatter plots are of the first two principal components plotted for haplotypic data, with the proportion of variance explained, from 6205 autosomal and 231 Z-linked loci (males only, because principal components analysis (PCA) does not accommodate heterogamy).

#### *2.2. ddRAD-seq Library Preparation*

Sample preparation for ddRAD sequencing followed the double-digest protocol outlined in DaCosta and Sorenson [54]. Genomic DNA (~1 μg) was digested with high fidelity versions of *Sbf*I and *Eco*RI restriction enzymes (New England Biolabs, Ipswich, MA, USA). Amplification and sequencing adapters containing unique barcode or index sequences were ligated to the sticky ends generated by the restriction enzymes. Libraries were size selected using gel electrophoresis (size range 300–450 bp) and purified using a MinElute Gel Extraction Kit (Qiagen) following the manufacturer's protocol. Size-selected fragments were amplified with Phusion high-fidelity DNA polymerase (Thermo Scientific, Pittsburgh, PA, USA) for 20 cycles, and purified using AMPure XP beads (Beckman Coulter, Inc., Indianapolis, IN, USA). Libraries were pooled in equimolar amounts determined via quantitative PCR (KAPA Biosystems, Wilmington, MA, USA). Single-end (150 bp) sequencing was completed on an Illumina HiSeq 4000 at the University of Oregon Core Genomics Facility. Raw Illumina reads are accessioned on National Center for Biotechnology Information (NCBI) Sequence Read Archive (BioProject PRJNA699594, Biosample accessions: SAMN17803887-SAMN17804091, see [55] for additional sample information).



#### *2.3. Bioinformatics*

Illumina reads were demultiplexed at the core facility and processed using the computational pipeline described by DaCosta and Sorenson ([54], custom Python scripts [56]). Briefly, the pipeline filters and clusters reads into putative loci based on sequence similarity (85%) using custom scripts and the UCLUST function in USEARCH v.5 [57]. Genomic positions of loci were determined by BLAST analysis [58] to the *Taeniopygia guttata* (Zebra Finch GenBank assembly reference GCA\_003957565.1) reference genome. MUSCLE v.3 [59] was used to align the reads within each cluster with genotyping of samples completed through custom scripts (See [60]). Genotypes were scored homozygous if > 93% of sequence reads were consistent with a single haplotype, whereas heterozygotes were scored if a second haplotype was represented by at least 29% of the reads, or if a second haplotype was represented by 20–29% of reads and the haplotype was present in other individuals. Loci were also "flagged" if the number of single-nucleotide polymorphisms (SNPs) was > 10, and if > 3 SNPs showed strong linkage. Alignments were visually inspected in Geneious (Biomatters Inc. San Francisco, CA), which allowed us to retain loci with insertion/deletions or high levels of polymorphism. To limit any biases due to sequencing error and/or allelic dropout, a minimum of 10 total reads was required to score a genotype as heterozygous with alleles with less than 5X coverage were scored as missing. Loci with a median depth of 10 per individual, < 10% missing genotypes, and < 10% flagged genotypes across all individuals were retained for downstream analyses.

Finally, autosomal and Z chromosome-linked loci were identified as described in Lavretsky et al. [60], with assignments based on differences in sequencing depth and homozygosity between males and females. Chromosomal positions across loci were attained by using blastn results against the *Taeniopygia guttata* reference genome from previous steps in the bioinformatic pipeline. These positions were used to verify marker type based on sequencing depth. Because females have only one Z chromosome, Z-linked markers in females are expected to appear homozygous and be recovered at about half the sequencing depth of males.

#### *2.4. mtDNA Sequencing*

Rusty blackbird individuals were sequenced at the mtDNA control region domain I and II. We amplified a 543-base pair (bp) fragment using primer pairs RUBL\_CR114L (5'–TCTTTGCCCCATCAGACAGC–3') and BBCR\_Rev1 [61]. Polymerase chain reaction (PCR) amplifications, cycle-sequencing protocols, and post-sequencing processing followed Sonsthagen et al. [62], with one exception: excess dNTPs and primer were removed using ExoSAP-IT (ThermoFisher Scientific, Waltham, MA, USA). PCR products were sequenced at Functional Biosciences, Inc. (Madison, WI, USA). For quality control purposes, we extracted, amplified, and sequenced 10% of the samples in duplicate. No inconsistencies in mtDNA sequences were observed between replicates. Sequence data are accessioned on GenBank (accession numbers: MW574845-MW574901; see [55] for additional sample information).

#### *2.5. Population Divergence and Nucleotide Diversity*

We calculated (1) nucleotide diversity (π) of each ddRAD locus (autosomal and Z-linked) and overall (all loci combined) and (2) composite pairwise estimates of relative divergence (*φST*) between sample locations using a custom Python script (out2phistA.py [63]).

Haplotype (*h*) and nucleotide (π) diversity were calculated for mtDNA in ARLEQUIN 2.0 [64]. Fu's *FS* [65] and Tajima's *D* [66] were calculated to test the hypothesis of selective neutrality and evidence of population fluctuations as implemented in ARLEQUIN. We applied critical significance values of 5%, which requires a *p*-value < 0.02 for Fu's *FS* [65]. An unrooted haplotype network for mtDNA loci was constructed in NETWORK 5.0.1.1 (Fluxus Technology, Suffolk, England 2019) using the reduced median method [67] to illustrate possible reticulations in the gene tree because of homoplasy or recombination. The degree of genetic divergence within rusty blackbirds was assessed by calculating overall and pairwise *FST* [frequency-based) and *ΦST* using a nucleotide substitution model [68] in ARLEQUIN.

#### *2.6. Population Structure—ddRAD*

Population structure was analyzed using four complementary methods with different underlying assumption requirements: (1) principal components analysis (PCA, nonparametric method) to identify major trends in the distribution of genetic variation; (2) maximum likelihood clustering analysis to estimate the number of underlying populations using individual SNPs in the program ADMIXTURE (parametric method); (3) fineRADstructure utilizing haplotypes (concatenation of all variable sites at each locus) to assess contemporary genetic relationships based on shared co-ancestry; and (4) estimate effective migration surfaces (EEMS, [69]) to identify regions that deviate from a null model of isolation-by-distance (IBD).

First, a PCA was implemented on Autosomal and Z-linked loci separately using haplotypic/allelic data and the dudi.pca function in the adegenet R package [70,71]. As PCAs require individuals to be either diploid or haploid, we only included males (in birds, the sex with two copies of Z chromosome) in the analysis of the Z chromosome loci. We plotted individuals relative to the first two principal components to determine the degree that genetically similar individuals cluster into distinct geographic groups.

Second, maximum likelihood estimates of population assignments across individuals were obtained with ADMIXTURE v.1.3 [72,73]. We used all autosomal bi-allelic SNPs with singletons (i.e., rare SNPs observed in only one individual) excluded and without a priori assignment of individuals to populations. First, SNPs were formatted for analyses using

plink [74], following steps outlined in Alexander et al. [75]. ADMIXTURE analysis was run with a 10-fold cross-validation, and a quasi-Newton algorithm employed to accelerate convergence [76]. Each analysis used a block relaxation algorithm for point estimation and terminated once the change (i.e., delta) in the log-likelihood of the point estimations increased by < 0.0001. To limit any possible stochastic effects from single analyses, we ran 100 iterations at each population of *K* (from *K* of 1–18). The optimum K was based on the average of CV-errors across the 100 analyses per K; however, additional K's were analyzed for further population structure resolution. We then used the program CLUMPP v.1.1 [77] to determine the robustness of the assignments of individuals to populations at each K. First, ADMIXTURE outputs were converted into CLUMPP input files at each K using the R program PopHelper [78]. In CLUMPP, we employed the Large Greedy algorithm and 1000 random permutations to estimate final admixture proportions for each K with per sample assignment probabilities (Q estimates, the log likelihood of group assignment) based on all 100 replicates per K.

Third, we used the fineRADstructure program [79] to cluster individuals into populations with indistinguishable genetic ancestry using a haplotype-based approach. FineR-ADstructure focuses on the most recent coalescent events (common ancestry) providing information on recent sample relatedness which can be informative regarding levels of contemporary gene flow. Samples were assigned to populations using 5,000,000 iterations sampled every 1000 steps with a burn-in of 500,000. We used 1,000,000 iterations of the tree-building algorithm to assess genetic relationships among clusters. Finally, the output was visualized using the R scripts, fineradstructureplot.r and finestructurelibrary.r [80].

We implemented the spatial method EEMS [69] to estimate effective gene flow (*m*) and genetic diversity (*q*) in order to identify areas across the breeding range that deviate from the null expectations of IBD. This method is based on a stepping-stone model where individuals are allowed to move between neighboring demes and gene flow rates can vary by locality. Expected genetic dissimilarity under the model depends on sample location and gene flow rates. Regions where genetic dissimilarity decays more quickly than expected are identified as barriers to gene flow or, conversely, corridors where genetic dissimilarity decayed more slowly than expected. A migration surface that correlates genetic variation with geography is interpolated to visualize potential barriers or corridors to movement. In addition, the model estimates an effective diversity parameter (q), which is the expected within-deme coalescent time and is proportional to average heterozygosity.

We used the same set of SNPs as in Admixture and calculated a dissimilarity matrix using bed2diffs R code included with the EEMS package [81]. An outer coordinate file was constructed using Google Maps API v3 [82] that included the species' entire breeding distribution [52]. Based on preliminary runs, we adjusted parameters, so the accepted proportion of proposals of variance was at least between 10% and 40%. We ran three independent analyses using 10,000,000 burn-in steps followed by 50,000,000 MCMC iterations sampled every 2000 steps for each deme size (100,250). We checked for convergence and visualized effective migration and diversity surfaces using the R package rEEMSplots [69].

#### *2.7. Effective Population Size—ddRAD*

Contemporary effective population size (Ne) was estimated from 6184 loci with NeEstimator v.2.1 [83] based on two methods: the linkage disequilibrium (LD) method [84] which tests for the nonrandom associations among alleles at different loci formed by genetic drift in small populations [85] and the molecular co-ancestry method [86] which evaluates the level of allele sharing among individuals. Further, we excluded rare alleles below a range of allele frequency values (Pcrit) from the linkage disequilibrium model to evaluate the effects of low-frequency alleles on Ne estimates. Variance in Ne estimates across a range of Pcrit values suggests a history of gene flow and/or the presence of first-generation dispersers, whereas stable Ne estimates are indicative of isolated populations [87]. We estimated Ne using a haplotype-based approach and Pcrit values between 0.01 and 0.09 and without a frequency restriction. Confidence limits (95%) were determined by jackknifing over loci. Ne was not estimated for sites with a sample size < 10.

#### *2.8. Hindcasted Paleo-Distributions of Breeding Rusty Blackbirds*

Finally, we evaluated the paleo-hindcasted maps of the potential LGM breeding range of rusty blackbirds by Stralberg et al. [53]. This was to assess whether population isolation during the LGM may have contributed to genomic structuring detected in our analyses of mtDNA and ddRAD. Models of rusty blackbird breeding density were developed by fitting boosted regression trees [88] to bioclimatic indices derived from current climate normals (1961–1990) as predictors of species abundance data derived from surveys conducted across the species' boreal range [89]. Climate variables were chosen from a set of bioclimatic indices [90] based on several criteria, including relevance to vegetation distributions, avoidance of extreme collinearity, and a preference for seasonal over annual variables when they showed high correlations. The final set of variables included extreme minimum temperature, chilling degree days below 0 ◦C, growing degree days above 5 ◦C, seasonal temperature difference, mean summer precipitation, climate moisture index [91], and summer climate moisture index.

Models were then fed inputs from downscaled paleo-climate projections for 21,000 years before present, according to two global climate models, Community Climate Model (CCM1) and Geophysical Fluid Dynamics Laboratory model (GFDL, [92]), to develop hindcast projections of the species' potential LGM breeding distribution [53]. We used projections from Stralberg et al. [53] to develop modified maps of potential LGM density that included areas thought to be covered by the Cordilleran and Laurentide ice sheets. This was to show where suitable habitats may have existed in unglaciated micro-refugia within the major ice sheets or along coastlines adjacent to known refugia now submerged under the sea (e.g., Bering Land Bridge, Grand and Georges Banks; [38]).

#### **3. Results**

#### *3.1. Bioinformatics—ddRAD*

We obtained over 294,699,715 million raw sequencing reads (median = 1,432,707 reads per individual, range 901,000–1,943,040) with a maximum 150 bp length. Initial exploration of genotyping results revealed that most loci were unambiguously genotyped across samples. We removed seven samples that were deemed to be of close familial relationship (e.g., siblings) based on preliminary PCAs and fineRADstructure results, and field notes (location, date, and age of individuals sampled). For the remaining 198 samples, a total of 6443 clusters (i.e., putative single-copy loci) met the depth/genotype threshold. Of these loci, 6436 passed automated checks for alignment quality or passed thresholds after manual edits that yielded 42,446 SNPs or insertion/deletion (polymorphic sites) from 6381 polymorphic loci. Of those, 6205 loci and 231 loci were assigned to autosomal and the Z chromosome, respectively. Final datasets comprised loci with a median sequencing depth of 118 reads per locus per individual (median range = 73−175 reads/locus/individual), and on average 98.5% (minimum of 80.0%) of alleles per individual per locus were scored.

#### *3.2. Population Divergence and Molecular Diversity*

Autosomal nucleotide diversity across the 6205 ddRAD loci was similar for all locations (0.0058–0.0070) with overall value of 0.0061 (Table 1). The highest percentage of loci with no variation (i.e., nucleotide diversity equals zero) was found within Yukon Territory, Canada (21.9%) and Vermont, USA (19.5%). Similar pattern was observed with Z-linked loci with overall nucleotide diversity of 0.0050 (range 0.0037–0.0049, Table 1). Among the 231 Z loci, Yukon Territory (39.4%) and Vermont (38.5%) had the highest percentage of non-variable loci. It should be noted that these populations also had the smallest sample size.

Overall, we uncovered relatively moderate levels of genetic differentiation across sample locations with Z-linked loci showing a 1.5× higher level of differentiation than autosomal (Figure 2; *φSTAutosomal* = 0.019; *φSTZ-linked*= 0.029). Autosomal and Z-linked loci *ΦST* values ranged from 0.002 to 0.082 and −0.016 to 0.144, respectively, with highest degrees of pairwise differentiation found between Newfoundland and other locations and between eastern and western locales (Figure 2). This pattern is reflective in the number of loci showing elevated pairwise divergence (*φST >* 0.1) which ranged from 0.2–30.3% of loci (12–1880 loci) across all comparisons (overall 0.8%, n = 52) for autosomal loci and from 0.9–33.3% for Z-linked loci (overall 5.2%, n = 12).

**Figure 2.** Pairwise *ΦST* (above diagonal) and *FST* (below diagonal) for mtDNA control region (**A**) and pairwise *ΦST* estimate for 6205 autosomal loci (above diagonal) and 231 Z loci (below diagonal (**B**). Darker colors indicate higher values. The northeastern region where higher levels of genetic differentiation were found is outlined.

We uncovered 56 unique mtDNA haplotypes characterized by 36 variable sites. Nucleotide and haplotype diversity were generally similar across sampled locations, except Cordova, Alaska which exhibited lower levels of diversity (n = 4, π = 0.0012, *h* = 0.333), and Nova Scotia/New Brunswick, Canada (n = 7, π = 0.0069, *h* = 1.000) and Manitoba (n = 7, π = 0.0063, *h* = 1.000) which exhibited higher diversity (Table 1). Tajima *D* was not significant, which is consistent with a hypothesis of selective neutrality of mtDNA. Fu's *Fs* was significantly negative for eight populations, suggestive of historical population growth.

Two main haplotype groups were observed in the mtDNA network, although separated by one or three variable sites depending on evolutionary pathway (e.g., presence of reticulations, Figure 3). The first group consisted of mainly Alaska locales with only 2 samples from eastern region and 17 haplotypes though only one haplotype was predominately represented. Conversely, the second group consisted of samples from both regions but was predominately comprised of central and eastern locations and 39 haplotypes with no one dominate haplotype. Genomic structure was uncovered with higher levels of differentiation estimated between Newfoundland and other sample locales, as well as between northeastern and Alaska locales. Similar results were observed with ddRAD loci; pairwise *ΦST* values ranged from −0.046 to 0.677 (overall *ΦST* = 0.147), and pairwise *FST* values ranged −0.014 to 0.438 (overall *FST* = 0.073, Figure 2).

**Figure 3.** Unrooted mitochondrial DNA haplotype median-joining network for rusty blackbirds. Size of circle is proportional to the frequency of each haplotype observed.

#### *3.3. Population Structure—ddRAD*

For PCA, plotting samples relative to the first two principal component axes based on autosomal loci recovered four main clusters that included (1) Newfoundland, (2) Nova Scotia and northeast United States, (3) Ontario, and (4) central Canada and Alaska (Figure 1). These four clusters were also retained but overlapped more when using only Z-linked loci. When PCA included samples only from Cluster 4 (central Canada through Alaska), Anchorage and Canada (Manitoba and Alberta) samples appear slightly differentiated from all other Alaska samples, although there is overlap in PC components (results not shown).

All possible K values were explored across ADMIXTURE analyses (Figure 4A). When K = 2, samples from northeastern locales (Clusters 1 and 2 in PCA) and Alaska and central Canada (PCA Cluster 4) were assigned with high probability to unique clusters. Ontario (PCA Cluster 3) was intermediate between (~70–80% assignment to Alaska cluster) the two main groups. When K = 3 or 4 are considered, the same overall pattern remained except central Canada (Manitoba and Alberta) and Ontario make up a third cluster with variable assignment probability (40–98%, see Supplementary Figure S1).

FineRADstructure revealed more sub-structuring than ADMIXTURE analyses, with individuals clustering mainly by geographic proximity (Figure 4B, Supplementary Figure S2). Locales in northeast North America had the highest shared co-ancestry values with samples being assigned to (1) Newfoundland, (2) Nova Scotia, and (3) northeast United States. Ontario samples were assigned to their own population but in agreement with ADMIXTURE, it shared higher co-ancestry with all other groups/populations indicating connectivity to western (overall higher with central Canada than Alaska) and northeastern locales. Samples from the central and western breeding range were primarily assigned to three main groups: (1) Alaska (excluding most Anchorage samples), (2) Anchorage, and (3) central Canada (Alberta and Manitoba). Unlike northeastern North America and Ontario where groups

were mutually exclusive, there was some admixture indicated by individuals being placed in a non-origin group within these central Canada and western populations, for example, Cordova, Alaska and Yukon Territory individuals being grouped with central Canada.

**Figure 4.** Population assignment analysis using the program (**A**) Admixture and (**B**) fineRADstructure. Average assignment probabilities for K = 2 inferred from bi-allelic single nucleotide polymorphisms in ADMIXTURE. FineRADStructure coancestry matrix indicating pairwise genetic similarity between individuals. Inferred populations are indicated by clustering in accompanying dendrogram and locations are indicated by color blocks for general region (solid blocks on right, color was chosen based on locale making up highest proportion of cluster) and by individual (single bars at left). Colors correspond to sampling locations indicated in Figure 1. Co-ancestry values were capped at 60 for illustrative purposes as only a few comparisons were above that value (see Supplementary Figure S2).

EEMS analysis highlighted regions of lower gene flow than expected under IBD. Migration surfaces were similar across deme sizes with one exception. Regions with reduced gene flow were (1) south-central Alaska, (2) Yukon Territory (only present with deme size 250), (3) Alberta, (4) eastern and western Ontario, and (5) Maritime provinces of Canada (Newfoundland, New Brunswick, and Nova Scotia; Figure 5). All these regions had high posterior probabilities (> 0.90) except for the western boundary of Ontario. These boundaries roughly correspond to the population clustering observed in fineRADstructure. High gene flow was characteristic across the remaining distribution.

#### *3.4. Effective Population Size—ddRAD*

Although Ne estimates were 10-fold lower for Anchorage, Alaska, and Alberta based on the linkage disequilibrium method, 95% confidence limits overlapped for all of the sampled sites as the upper bounds were infinity (Figure 6). Variation in Ne estimates across Pcrit values was observed for Alaska sites (Tetlin and Yukon Flats), Ontario, and New Hampshire, indicative of past gene flow or the presence of first-generation dispersers affecting Ne estimates. Point estimates for Ne for Newfoundland were infinity, indicating there is no evidence that the population is not very large. However, lower bounds of 95% confidence levels using jackknife method ranged from 83.5 (Pcrit = 0.09) to 372.9 (Pcrit = 0.01) providing a plausible limit for Ne [93]. Conversely, Ne estimates based on the molecular co-ancestry method were lower for two Alaska sites (Anchorage and Tetlin) and Ontario (Table 1).

**Figure 5.** The estimated effective migration surfaces between all sampling locales (black circles) for deme size 250. White areas indicate gene flow rates consistent with isolation by distance expectations whereas shaded areas have dispersal rates that are higher (blue, corridors) or lower (orange, barriers) than expected under isolation-by-distance (IBD). We note that orange area around the Yukon Territory sampling location showed opposite pattern when deme size was lower than 250 (higher than average gene flow rate). For all other areas deme size did not change results.

**Figure 6.** Effective population size (Ne) estimates as a function of excluding rare alleles (Pcrit) in rusty blackbirds sampled across their North American distribution. Point estimates of Ne are log transformed with values >5 signifying infinitely large estimates. Colors correspond to Figure 1.

#### *3.5. Hindcasted Paleo-Distributions of Breeding Rusty Blackbirds*

Suitable climate conditions for breeding rusty blackbirds hindcasted to the LGM by Stralberg et al. [53] were located primarily in what is now the central and eastern United States, with some potential suitable habitat in northwestern regions of the United States and Canada (Figure 7). Projected LGM densities were lower, and distributions were more limited and fragmented compared to the species' current distribution and abundance patterns. Despite substantial variation between them, both CCM1 and GFDL models hindcasted suitable climates for paleo-populations of rusty blackbirds in or adjacent to glacial refugia in (1) Alaska and Yukon Territory (Beringia), (2) insular British Columbia (Vancouver Island and Haida Gwaii), (3) Newfoundland (Grand Banks), and (4) New England and the Canadian Maritimes (Georges Banks), as well as south of the ice sheets in (5) British Columbia and Washington and (6) the Great Plains, Upper Midwest, and Northeast United States (Figure 7).

**Figure 7.** Predictions of current breeding density (**left**) versus hindcasted paleo-breeding densities of rusty blackbirds (based on [53]. Paleo-breeding densities were hindcasted by fitting models of current rusty blackbird density with bioclimatic indices downscaled by Roberts and Hamann [92] for the last glacial maximum (21,000 YBP) using two global climate models: Community Climate Model (CCM1, **middle**) and Geophysical Fluid Dynamics Laboratory (GFDL, **right**). Model-predicted density values range from 0 to 1.4 pairs/ha. Overlaid on hindcast projections are level 1 ecoregion boundaries in black [94] and the extent of Last Glacial Maximum (LGM) ice sheets in transparent gray [95,96].

#### **4. Discussion**

The geographic patterns in genomic structure we detected across the rusty blackbird's breeding range conformed to our hypotheses. (1) An east–west partition in genomic structure was observed between rusty blackbirds nesting in the western and central boreal regions versus the eastern boreal region. This was consistent with the east–west migratory divide detected for the species using stable isotopes [35] and suggests long-term separation of populations. (2) As expected, based on insular isolation and plumage differences [37], birds in Newfoundland were differentiated from birds from other sampled sites. (3) Populations within both eastern and western regions exhibited subtle genomic structuring and restricted gene flow, indicating dispersal is limited by discontinuities in habitat, physical barriers, philopatry, or migratory behavior. Further, Ontario appears to be an area of secondary contact between birds originating from eastern and western lineages identified in ADMIXTURE and fineRADstructure analyses. Together, these results indicate that historical and contemporary processes are shaping the distribution of genomic variation among populations of rusty blackbirds across their boreal distribution.

#### *4.1. Pleistocene Influences on Patterns of Genomic Diversity*

While the species' current nesting distribution is largely contiguous across the boreal forest biome, the distribution hindcasted to 21,000 years before present was displaced

south (other than portions of Beringian Alaska) and was fairly discontinuous and limited (Figure 7, [53]). During the LGM, the glacial ice sheets covered most of northern North America, except Alaska and areas along Pacific and Atlantic coasts [97], and may have sundered the rusty blackbird's nesting range, isolating populations in separate western and eastern refugia, and promoting the partition in genomic variation we detected in multiple analyses of mtDNA and ddRAD (autosomal and Z-linked) loci. In the same way, glacial vicariance is hypothesized to have led to Pleistocene speciation in several sister pairs of temperate conifer boreal bird species [47,98,99]. In addition, several passerines with trans-boreal distributions share this east–west divide in genomic diversity with secondary contact between divergent lineages occurring in the central boreal: mixing occurs from Alberta to Manitoba within the yellow warbler (*Setophaga petechia*, [100,101]), from Alberta to Ontario within Wilson's warbler (*Cardellina pusilla*, [15,102]), from Saskatchewan to Manitoba within the Canada jay (*Perisoreus canadensis*, [41,44]), from Alberta to Ontario within the golden-crowned kinglet (*Regulus satrapa*, [103]), and from Manitoba to Ontario within the rusty blackbird (this study). However, the east–west divide within the blackpoll warbler (*Setophaga striata*) was attributed to isolation by distance and not glacial vicariance [46]. These concordant breaks in genomic diversity across multiple trans-boreal species emphasize the strong influence that the Pleistocene ice sheets played in shaping how genomic variation is arrayed across northern North America in boreal avifauna.

Models of the paleo-breeding distribution indicated that the nesting habitat for rusty blackbirds could have been present in four potential glacial refugia: (1) Alaska (Beringia), (2) Atlantic Shelf, and south of the ice sheets in (3) western (Cordilleran) United States, and (4) eastern (Laurentide) United States, with the eastern region likely a core area based on the relatively high densities inferred during the LGM [53]. These four regions coincide with the locations of glacial refugia proposed for the boreal chickadee [42], the Canada jay [41,44], and the black spruce (*Picea mariana*, [104])—the tree species most often selected for nesting by rusty blackbirds [105]. Rusty blackbirds and other spruce-associated species may have therefore followed the post-glacial colonization routes inferred for black spruce from pollen, fossils, genetics, and ecological modeling [41,106,107]. The spatial apportionment of genomic diversity does suggest that rusty blackbirds occupied at least two refugia during the LGM, although it is not clear from model hindcasts which regions were the refugia nor from the genetic data which sample locations represent refugial populations. Recent colonization of deglaciated areas, whether via long-distance dispersal or leading-edge expansion from glacial refugia, leaves predictable signatures, notably reduced genomic diversity associated with founder events followed by founders preventing subsequent waves of colonizers [108]. Levels of genomic diversity, however, were similar across the rusty blackbird distribution (Table 1). Dispersal likely continued from refugial populations into founding populations in deglaciated areas either via short movements or continued long-distance dispersal. Connectivity between refugial and founding populations would have maintained genomic diversity because effective population sizes would not have been markedly reduced [39,109], thereby erasing the genetic legacy of founder events associated with post-glacial colonization. Further, the eastern clade had high haplotype diversity with no single dominate haplotype suggesting that Newfoundland, Canada Maritime provinces, and New England states were colonized by rusty blackbirds originating from the "core" eastern refugium and possibly the Atlantic Shelf refugium as the effective population size would need to be large to retain and maintain genetic diversity through the LGM. Conversely, the mtDNA network for the western clade had a star-like pattern with a single dominant haplotype predominately represented by Alaska birds. This suggests that the western breeding range from Alaska to Manitoba was colonized by a small refugial population of rusty blackbirds expanding their range from a western or Alaska (Beringia) refugium that supported lower nesting density, and therefore presumably lower effective population size during the LGM.

#### *4.2. Isolation in Newfoundland*

Rusty blackbirds occupying Newfoundland were genetically divergent from all sampled locations across marker types. The presence of genetic structure is concordant with the putative subspecies, *Euphagus carolinus nigrans*, described by Burleigh and Peters [37] to nest on Newfoundland and Magdelen Islands, winter in South Carolina, and have a distinctive plumage that is darker, glossier, and bluer than the nominate form. Burleigh and Peters [37] described an additional 12 endemic subspecies of passerines breeding on Newfoundland that also had darker plumages than their mainland counterparts. Similar to rusty blackbird, several of these species and others have since been found to have genetically distinct populations on Newfoundland: the American redstart (*Setophaga ruticilla*, [45]), the blackpoll warbler [46], the boreal chickadee [42], the purple finch (*Haemorhous purpureus*, [110], the Canada jay [44], the gray-cheeked thrush [111], and the black-capped chickadee (*Poecile atricapillus*, [112]). Cabot Strait (≥104-km wide) and the Strait of Belle Isle (≥15-km wide) separating Newfoundland from mainland Canada therefore seemingly act as strong physical barriers to dispersal by the rusty blackbird and several other passerines.

While Newfoundland rusty blackbirds were genetically differentiated, they were peripheral on the mtDNA network and shared several mtDNA haplotypes with birds breeding elsewhere in eastern North America. This coupled with the limited structure we observed in ADMIXTURE and other analyses of ddRAD loci indicate that divergence between Newfoundland and other eastern locales likely arose post-Pleistocene and has been maintained through restricted dispersal. The observation of shallow divergence is consistent with several other passerines that nest on Newfoundland [44,46,110]. In contrast, other species show much deeper genetic divisions between birds in nesting areas along the northern Atlantic coast versus nearby locales [40,42,112]. In these species, it is more plausible that Newfoundland populations were isolated during the Pleistocene on the Atlantic Shelf refugium offshore of Newfoundland along the Grand Banks [38,97]. Although hindcasting models indicate that suitable habitat was available in Newfoundland during the LGM, the Laurentide ice sheet covered most of the region. Therefore, if rusty blackbirds currently occupying Newfoundland were present in the Atlantic Shelf refugium (Grand Banks, [38,97]), individuals were likely restricted into small isolated area(s) with low numbers. As genetic drift is a strong force shaping genomic diversity in small isolated populations, the lack of deep partitions in genomic variation suggests rusty blackbirds colonized Newfoundland post-Pleistocene and that genomic and morphological variation likely evolved recently as shown for the widely distributed and phenotypically diverse *Junco* species complex [113].

#### *4.3. Contemporary Influences on Patterns of Genomic Variation*

We also found several lines of evidence that contemporary processes are limiting dispersal of rusty blackbird populations. First, the maintenance of two distinct genetic groups in western versus eastern North America is indicative of continued restrictions to dispersal between regional nesting areas. This genetic divide mirrored the general migratory divide between western and eastern flyways identified with stable isotopes [35]. Thus, differences in migratory pathways or timing of migration may reinforce reproductive isolation or spatial segregation of regional rusty blackbird populations in the same way as suggested for other passerines with northern distributions [114–117]. We also identified a potential area of secondary contact between eastern and western lineages of rusty blackbirds from a single sample location at the northern border between Ontario and Quebec which had similar levels of recently shared ancestry with both western and eastern regions. Genetic samples coupled with tracking studies from additional eastern locales would help determine the geographic extent of the mixing zone and the degree that migratory behavior may be restricting genomic connectivity [115].

Second, we found evidence of restricted dispersal within eastern and western North America in the form of subtle structure in ddRAD loci among most sampled locales (e.g., Nova Scotia versus other northeastern locales and Alberta/Manitoba versus Alaska) that

deviated from expectations based on IBD on a regional scale (Figure 5). Furthermore, three areas (Anchorage, Alberta, and Newfoundland) had similar estimates of effective population size (Ne) across a range of rare alleles, suggesting genomic diversity in these populations is not influenced by ongoing gene flow from adjacent populations. While mountain ranges and ocean bodies are obvious barriers isolating rusty blackbirds in Anchorage and Newfoundland, respectively, there are no clear physical barriers restricting dispersal between Alberta and Manitoba within central Canada. Despite the lack of physical barriers, Alberta individuals were assigned to a non-origin grouping (24%, 5/21) less often than Anchorage individuals (38%, 9/24) in the fineRADstructure analysis, suggesting that other factors are impeding dispersal. Indeed, rusty blackbird nest in boreal forest wetlands that are often distributed in discrete patches separated by unsuitable upland habitat or fragmented by frequent natural disturbance [104,118]. Although rusty blackbirds as well as other passerines migrate over long distances, many species return to near their natal sites to nest where local landscape features influence movements to locate new nesting areas [119]. Other forest dependent birds have shown a reluctance to disperse across large gaps of non-forested habitat to locate new nesting grounds [14,51,120,121]. This type of behavior, if exhibited by rusty blackbirds, may limit dispersal and contribute to genomic structuring among some rusty blackbird nesting areas.

#### **5. Conclusions**

Across the breeding range, rusty blackbirds exhibited genomic structuring evident of restricted gene flow, which may limit the species' adaptive capacity to respond to rapid environmental change. The North American boreal biome is a mosaic of wetland complexes and forests that are projected to be transformed as the Earth's climate continues to warm and increase the frequency and magnitude of boreal disturbances such as drought, permafrost thaw, fire, and insect outbreaks [122]. Under various simulations of climate-mediated ecological change over the 21st century, the boreal biome is projected to contract by up to 42% [123], and boreal birds are projected to both dramatically shift their ranges northwards and upwards in elevation and suffer disproportionately high losses in population size and range extent among North America avifauna [22,89,122,124,125]. The rusty blackbird is particularly vulnerable to projected reductions in suitable breeding habitat, which could result in the loss of more than half of the species' breeding range [125] and population numbers [89]. These future declines will exacerbate the species' already steep global population decline and southern range retraction since the mid-20th century [25,26]—the latter already linked to regional trends in warming [33]. Additional research on genotype-environment associations using functionally relevant loci (e.g., transcriptome or gene expression analysis) can build off of the foundation of this study to identify breeding areas that may be more vulnerable to stochastic events as well as areas that pose high conservation value for the species as the climate continues to change (e.g., [101,126]).

The rusty blackbird has an immense migratory range (breeding across the continental boreal biome, wintering over the eastern half of United States) and the many stressors suspected to be contributing to its decline are hypothesized to vary widely across breeding areas and over the annual cycle [26,27,127]. Efforts to understand the causes of decline and efficiently link conservation across this species' annual cycle will therefore benefit from a more comprehensive knowledge of migratory connectivity than the general east–west migratory divide identified through stable isotope and tracking studies [34–36]. Our study is a foundational step in gaining this knowledge as it provides a basis for researchers to infer the natal origins of birds sampled at key migration stopover sites and important wintering areas (e.g., [15]). Understanding migratory connectivity across the rusty blackbird's nonbreeding range would, for example, allow researchers to weigh the relative contributions of summer versus winter environmental change on vital rates and population trends (e.g., [128,129]) and enable wildlife managers to strategically target habitat restorations throughout the annual cycle for genetically distinct populations [130]. As the boreal avifauna is among the most rapidly declining groups of birds in North America [23], the

integration of information on connectivity with the science and management of recovering rusty blackbird populations could serve as a model for how to restore other poorly studied declining boreal species [27].

**Supplementary Materials:** The following are available online at https://www.mdpi.com/1424-281 8/13/3/103/s1, Figure S1: Admixture results for population clustering, Figure S2: FineRADStructure co-ancestry matrix.

**Author Contributions:** Conceptualization, R.E.W., S.A.S., L.L.P., S.M.M., J.A.J., D.W.D.; resources, S.M.M., J.A.J., L.L.P., D.W.D.; methodology, R.E.W., S.A.S.; formal analysis, R.E.W., S.A.S., D.S.; investigation, R.E.W., S.A.S.; writing—original draft, R.E.W., S.A.S.; writing—review and editing, R.E.W., S.A.S., L.L.P., S.M.M., J.A.J., D.S., D.W.D.; funding Acquisition, J.A.J., D.W.D., S.A.S.; project administration, S.A.S., L.L.P., S.M.M., J.A.J., D.W.D.; data curation, R.E.W., S.A.S.; fieldwork, J.A.J., L.L.P., S.M.M., All authors have read and agreed to the submitted version of manuscript.

**Funding:** This research was funded by the U.S. Fish and Wildlife Service, Division of Migratory Bird Management and U.S. Geological Survey, Wildlife Program of the Ecosystems Mission Area.

**Institutional Review Board Statement:** Sample collections were approved by Institutional Animal Care and Use Committees of the U.S. Fish and Wildlife Service, Alaska Department of Fish and Game (Alaska: IACUC #2008004, approved May 2008) and University of Alberta (Alberta: #AUP00001523).

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data used in the present study are accessioned in Sonsthagen et al. [55], NCBI Sequence Read Archive (BioProject PRJNA699594, Biosample accessions: SAMN17803887-SAMN17804091) and NCBI nucleotide database (MW574845-MW574901).

**Acknowledgments:** We thank Erin Bayne, Paulson Des Brisay, Matthew Brooks, Terry Chesser, Andrew Curtis, Niels Dau, Lucas DeCicco, David Evers, Carol Foss, April Harding-Scurr, David Loomis, Bruce Rodrigues, Pam Sinclair, David Shaw, Marian Snively, David Tessler, and James Wright for their critical contributions to sample collection. Sample collections in Alaska were approved by Institutional Animal Care and Use Committees of the U.S. Fish and Wildlife Service and Alaska Department of Fish and Game. This research used resources provided by the Core Science Analytics, Synthesis, and Libraries (CSASL) Advanced Research Computing (ARC) group at the U.S. Geological Survey. Any use of trade, product or firm names is for descriptive purposes only and does not imply endorsement by the US Government. The findings and conclusions in this article are those of the authors and do not necessarily represent the views of the USFWS. This article has been peer reviewed and approved for publication consistent with USGS Fundamental Science Practices.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


*Article*
