Habitat‐linked genetic structure for white‐crowned sparrow (Zonotrichia leucophrys): Local factors shape population genetic structure

Abstract Ecological, environmental, and geographic factors all influence genetic structure. Species with broad distributions are ideal systems because they cover a range of ecological and environmental conditions allowing us to test which components predict genetic structure. This study presents a novel, broad geographic approach using molecular markers, morphology, and habitat modeling to investigate rangewide and local barriers causing contemporary genetic differentiation within the geographical range of three white‐crowned sparrow (Zonotrichia leucophrys) subspecies: Z. l. gambelii, Z. l. oriantha, and Z. l. pugetensis. Three types of genetic markers showed geographic distance between sampling sites, elevation, and ecosystem type are key factors contributing to population genetic structure. Microsatellite markers revealed white‐crowned sparrows do not group by subspecies, but instead indicated four groupings at a rangewide scale and two groupings based on coniferous and deciduous ecosystems at a local scale. Our analyses of morphological variation also revealed habitat differences; sparrows from deciduous ecosystems are larger than individuals from coniferous ecosystems based on principal component analyses. Habitat modeling showed isolation by distance was prevalent in describing genetic structure, but isolation by resistance also had a small but significant influence. Not only do these findings have implications concerning the accuracy of subspecies delineations, they also highlight the critical role of local factors such as habitat in shaping contemporary population genetic structure of species with high dispersal ability.


| INTRODUC TI ON
Dispersal of individuals from one habitat patch to another is a simple concept, yet dispersal can have a profound effect on gene flow, genetic diversity, adaptation, and population dynamics (Bowler & Benton, 2005;Olah et al., 2017;Ronce, 2007). This is especially true for species that exhibit habitat specialization or low mobility where environmental and landscape features, including climate, habitat, and elevation, can act as barriers to dispersal (Geffen et al., 2004;Kershenbaum et al., 2014). For these reasons, the influence of landscape on population connectivity and genetic structure is complex and depends on individual behavior and the landscape composition (Coulon et al., 2004;Holderegger & Wagner, 2008;Lee-Yaw et al., 2009;Rissler, 2016).
The relationship between gene flow and dispersal potential is not always apparent. Many species capable of flight are thought to have high population connectivity simply because they can cross or circumvent physical barriers such as anthropogenic development, large bodies of water, or mountain ranges, but this is not always true (Holderegger & Wagner, 2008). With the advent of landscape genetics, more studies show genetic differentiation can occur within contiguous populations due to microclimates and behavioral differences unrelated to physical barriers (Dubay & Witt, 2014;Engler et al., 2014;Porlier et al., 2012). For example, mitochondrial DNA differentiation exists between rufous-collared sparrow (Zonotrichia capensis) populations locally adapted to low versus high elevation microclimates on the same mountain slope (Cheviron & Brumfield, 2009). Similarly, in the absence of bluegill sunfish (Lepomis macrochirus), some pumpkinseed sunfish (Lepomis gibbosus) changed from preying on snails to feeding on Daphnia, which resulted in morphological differences between the ecotypes (Robinson et al., 1993). Therefore, even in the presence of high dispersal potential and in the absence of conspicuous physical barriers, population genetic structure can occur at small spatial and short temporal scales (Porlier et al., 2012).
Physical barriers such as mountains are an effective barrier to dispersal for many species (Fjeldsâ et al., 2011;Hooge, 2003;Robertson et al., 2009). However, if a species with breeding sites on either side of a mountain range share a wintering ground and have low philopatry, this could lead to the incorrect conclusion that the mountains do not act as a barrier. In addition to physical barriers, environmental and habitat variables also affect genetic variation F I G U R E 1 Breeding season range map of white-crowned sparrow subspecies, Z. l. gambelii, Z. l. leucophrys, and Z. l. oriantha in shades of blue, Z. l. pugetensis in red, and Z. l. nuttalli in orange. Areas with cross-hatching indicate areas of overlap between subspecies. Adapted from Dunn et al. (1995). Example photographs of subspecies show the following: (a) Z. l. gambelii has pale lores (small feathers between bill and eye) and orange bill. (b) Z. l. leucophrys has a thin line of black in the lores near the eye and a light pink bill. (c) Z. l. oriantha is very similar in appearance to Z. l. leucophrys, but has flared crown stripes and heavy black lores, a dusky pink bill with darker red on culmen (upper bill) and on average has a smaller body mass, longer wings, and shorter tarsi. (d) Z. l. pugetensis has pale lores and a yellow bill. Z. l. nuttalli is virtually indistinguishable in the field from Z. l. pugetensis (Dunn et al., 1995;Morton, 2002). Photo credits to C. Welke © and Don Robinson © (a, b, and d; with permission) TA B L E 1 Number of samples sequenced (n) at each sampling site and for each subspecies and location (ID) for control region (CR) and Aldolase B (AldoB6) sequences and microsatellite genotypes     (Sexton et al., 2014). For example climate, habitat and diet have been found to influence genetic variation for a wide variety of aquatic and terrestrial species (Pilot et al., 2006;Selkoe et al., 2010;MacDonald et al., 2020). When gene flow is restricted between habitats or ecosystems, it can lead to the formation of genetically distinct ecotypes (Kess et al., 2018;Parchman et al., 2016).
The white-crowned sparrow (Zonotrichia leucophrys) is an ideal species for testing predictions about the effects of rangewide and local barriers on genetic structure. It is a widespread North American passerine whose range spans multiple barrier types and is composed of five subspecies with different song dialects, habitat preferences, phenotype, and migration behaviors (Figure 1).

F I G U R E 2
Map of the 15 sampling sites in this study and minimum spanning networks for CR using sequences from all sites except CNP and 66 birds (a) and AldoB6 using sequences from 12 sites and 54 birds (85 sequences after phasing male haplotypes) (b). Each circle represents a haplotype, with the number of dashes across connecting lines showing number of base pair differences. Colors correspond to the haplotype's population of origin. A 9 bp insertion linked to a 19 bp deletion in some AldoB6 haplotypes is indicated by "indel." The single sample from OK is displayed in the network but omitted from F ST comparisons. Site abbreviations available in Table 1 F I G U R E 3 STRUCTURE plot of all populations with optimal number of genetic clusters of K = 4 as determined by highest log-probability and ΔK for 328 samples (a). (b) Hierarchical analysis of the first group (yellow) had an optimal number of K = 2. (c) Substructure according to ecosite type for northern samples also had K = 2. Bayesian analysis using TESS v2.3 yielded an optimal K = 3 in the overall clustering (d), K = 2 for each of cluster 1 (light green), cluster 2 (light blue), Colorado elevation groups (purple) (e), and northern ecosite groups (f). Above each plot is a bar showing subspecies of each individual as Z. l. gambelii (gray), Z. l. oriantha (black), Z. l. pugetensis (striped), or hybrid gambeliioriantha (yellow). Unknown subspecies are uncolored Previous studies suggest levels of gene flow differ among subspecies (Austen & Handford, 1991;Chilton et al., 1990;Lein & Corbin, 1990;Soha et al., 2004) and hybridization has been reported between two Z. l. gambelii and Z. l. oriantha (Lein & Corbin, 1990). To date, no studies have extensively examined the factors contributing to genetic variation within this species. Do genetic patterns arise as a result of environment or as a result of isolation by distance (Sexton et al., 2014)? In this study, we use genetic and environmental data for white-crowned sparrows to examine how landscape composition and habitat mediate gene flow.
We examined population genetic structure for three subspecies Z. l. pugetensis, Z. l. oriantha, and Z. l. gambelii of the white-crowned sparrow. We sampled the northern Rocky Mountains, a broad north-south belt of mountains interspersed with rift valleys and glacier-carved basins (Cruden & Hu, 1999;Holland, 1976 Rivers et al., 2019).

| Genetic diversity analyses
Out of 54 birds used for AldoB6, 31 were male. Haplotypes from these males were reconstructed using PHASE in DnaSP v5 (Rozas et al., 2017) (Pritchard et al., 2000). The program was run with correlated allele frequencies and the admixture model with sampling locations as locpriors. Ten independent runs were performed with 50,000 burn-ins and 200,000 McMC repetitions for each K value from K = 1-10. To determine the optimal number of clusters, STRUCTURE HARVESTER v0.694 (Earl & vonHoldt, 2012) was used to average values of LnPr(X|K) and ΔK for each K (Evanno et al., 2005). Seven of the ten northern populations showing admixture were run separately with the same settings to determine whether additional structure was present. To determine whether genetic structure corresponded to ecosite type (northern populations) or elevation (CO and OR), we used ecosite or elevation as locpriors and ran STRUCTURE using the same settings as above.
To complement our STRUCTURE analysis, we used a second Bayesian clustering algorithm TESS v2.3 (Chen et al., 2007) that incorporated spatial information. TESS was run for K values from 1 to 10 using 5,000 burn-ins and 100,000 sweeps, with a spatial interaction parameter (Ψ) of 0.6. The optimal K was determined by choosing the K where the deviance information criterion (DIC) value began to plateau (Chen et al., 2007). As with STRUCTURE, the same groupings for investigating admixture, ecosite types, and elevation were run separately for K values from 1 to 6 with the same settings to uncover possible substructure.

| Multivariate analyses
To further assess genetic structure from a multivariate perspective,

Principal Coordinate Analysis (PCoA) and a Principal Component
Analysis (PCA) were performed in R (R Core Team 2021). For the PCoA, a matrix of multivariate genetic distances (F′ ST values) between all populations was plotted against their geographic coordinates. A three-dimensional PCoA graph was made using the Scatterplot3D package (Ligges & Mächler, 2003) to visualize the first three principal coordinates. In addition to comparing genetic differences, we also used multivariate analyses to examine morphological differences among ecotypes. We conducted a PCA on six morphological measurements: bill length, bill depth, bill width, tarsus length, uncompressed wing chord, and body mass. For this analysis, we compared morphological differences from three ecotypes: alpine coniferous, riparian deciduous, and disturbed-gas site. Morphological differences were summarized with a standard multivariate ordination using the ggbiplot package (Wickham, 2009).

| Species distribution models
We constructed a species distribution model (SDM) to display potential dispersal routes of white-crowned sparrows using least-cost path (LCP) and least-cost corridor (LCC) models. The SDM was made in ArcMap 10.1 (ESRI, Redlands, CA) using SDMtoolbox v1.1 (Brown, 2014

| Correlates of genetic structure with distance and dispersal
To model the most likely dispersal routes and dispersal costs, LCP and LCC analyses were conducted with SDMtoolbox v1.1c (Brown, 2014

| Isolation by distance and resistance
We used Mantel and partial-Mantel tests to examine the effects of isolation by distance (IBD) and isolation by resistance (IBR) on genetic differentiation. All Mantel and partial-Mantel tests were conducted in GenoDive 3.04 (Meirmans & Van Tienderen, 2004). We converted pairwise F′ ST values to a distance measurement using the formula F′ ST /(1 − F′ ST ) and calculated geographic distances between sites in GenAlEx 6.05 (Peakall & Smouse, 2012). For IBR, we generated a least-cost path distance matrix using the along-path cost of the LCP calculated in ArcMap 10.1. The along-path cost is the total sum of the friction values that characterize the LCP and allows for the testing of IBR between sample sites (Etherington, 2011

| Genetic diversity
The mitochondrial CR had a total of 17 haplotypes; six haplotypes were shared among multiple populations, while the remaining 11 haplotypes were restricted to a single population ( differentiation between all northern population pairwise comparisons, as well as southern population pairwise comparisons (range: 0-0.30). All three subspecies were genetically distinct based on control region sequences (range: 0.15-0.40). We observed low to high genetic variation for comparisons of AldoB8 sequences (range: 0-0.72). For the 15 AldoB8 θ ST comparisons within northern populations, 46% were significant, and 20% were significant within the 10 southern population comparisons (Table 2). In the north, BV and MK accounted for all the significant differences and OR-L in the south.
Similarly, between the north and south, MK, LE, and OR-L accounted for all but one of the significant results (Table 2). Among the three subspecies, Z. l. oriantha and Z. l. gambelii were significantly different from Z. l. pugetensis, but not from each other. Genetic divergence at this marker was relatively low (range: 0.04-0.12).
We observed high population genetic structure based on micro-     The PCA using morphometric data from all northern populations showed two distinct ecosite groupings similar to STRUCTURE and TESS with the first and second axes accounting for 32.8% and 18.7% of the variation, respectively. When grouping individuals by ecosites, the six phenotypic traits clustered into AC and RD/D-G ecosites.

| Population genetic structure
The samples from the disturbed group only include the D-G population in CNP (Figure 4).

| Species distribution models and potential dispersal routes
The contemporary SDM for each subspecies closely followed their known distributions in North America (Figures 1 and 5)

| Correlates of genetic differentiation with distance and dispersal
For our models examining all populations excluding the two Oregon populations, IBD and IBR were significant (range: r = .24-.29, p ≤ .05).
We were unable to separate the effects of these two processes, IBD

| D ISCUSS I ON
We compared the influences of barriers and ecotype on population genetic structure at various spatial scales, and we hypothesized that geographic distance and mountains would act as barriers to gene flow within and among three subspecies of white-crowned sparrows.
Although IBD and IBR influence rangewide patterns, genetic patterns

F I G U R E 4 A standard principal component analysis of 113
individuals with multivariate analysis of six morphological features (wing length; bill width, length, and depth; tarsus length; and body weight). Ellipses show individuals from alpine coniferous (light green), riparian deciduous (dark green), and disturbed habitat (blue) ecosites as compared in TESS and STRUCTURE (Figure 3). Disturbed habitat individuals are from the disturbed-gas plant site only, as no morphological data were available for the 11 townsite individuals, and CO and OR populations were omitted at local scales were not influenced by these factors. Our results indicate a role for local ecological factors and microclimate acting as barriers and influencing both rangewide and local genetic patterns.

| Macrogeographic barriers: distance and hybrid zones
Physical distance played a large role in describing the rangewide genetic differentiation in Z. l. gambelii and Z. l. oriantha (Figure 2, Table 4). Within each of these ranges, large amounts of suitable habitat are present; however, habitat and dispersal are not homogenous as evident by the genetic analyses. Both nuclear datasets showed north-south differences and evidence of dispersal across mountain ranges (Table 3,  F I G U R E 5 Contemporary SDM created using occurrence data from three subspecies of white-crowned sparrow, environmental variables, a vegetation cover layer, and elevation layer. Areas of the most suitable environmental and habitat conditions (i.e., ecological niche) for each subspecies are in warmer colors (left). Least-cost corridor (LCC) projections of dispersal routes between sites within the breeding ranges of Z. l. gambelii and Z. l. oriantha based on the SDMs with base maps showing higher elevations in lighter gray. Dispersal costs are coded by red representing areas with low resistance (i.e., low dispersal cost), and blue representing high resistance (right). The models were created using the SDM toolbox (Brown, 2014) in Arc GIS® and MaxEnt (Phillips et al., 2006) TA B L E 4 Mantel tests comparing genetic distance with both geographic distance and dispersal resistance Note: F′ ST values are compared against Euclidian distances between populations for a test of isolation by distance (IBD), and against the least-cost path resistance values for a test of isolation by resistance (IBR). Significant p-values are in bold. Variables following vertical lines indicate the variable that was controlled for during partial-Mantel tests. Hooge, 2003). Similarly, the genetic similarity of BV and RV could potentially be explained by dispersal through the Rocky Mountain Trench.
The differentiation of our CO samples from the northern individuals using microsatellite data may be due to the sampling gap or hybrid zones. Within the Rocky Mountain area of AB and BC, a hybrid zone exists between Z. l. oriantha and Z. l. gambelii, whereas in CO they are pure Z. l. oriantha (Figure 1). In the sampling gap between northern and CO sparrows, there could be patterns of introgression defining the hybrid zone edge. Alternatively, if there is clinal variation in ecosystems, the differentiation of northern populations of CO could be an artifact of that sampling gap (Barton & Hewitt, 1989). Exploring these two possibilities is important because the patterns of population differentiation within the hybrid zone do not correspond to subspecies, but to ecosystems.

| Microgeographic barriers: Forest type and disturbed ecosystems
We found genetic clusters corresponding to alpine coniferous (AC) and riparian deciduous (RD) ecosite types and individuals from disturbed ecosites grouped with the riparian deciduous individuals ( Figure 3). Morphological and genetic differences have been observed across heterogeneous landscapes in plants (Gram & Sork, 2001), insects (Hamer et al., 2003), aquatic invertebrates (Etter et al., 2005), and terrestrial vertebrates (Barley et al., 2015). Local environmental conditions and ecosystem have an important role in determining population genetic structure of white-crowned sparrows. Individual white-crowned sparrows from the same sampling area in OK belong to different genetic groups corresponding to the ecosite they were sampled in, although due to the low number of individuals collected from this area, our results should be viewed conservatively. No reports exist on the contemporary status of white-crowned sparrows in the OK area. Brooks (1909) and Kermode (1904) report Z. l. gambelii breeding in this area, but other studies (Krannitz, 2007;Krannitz & Rohner, 2000) claim white-crowned sparrows are transient. Our museum samples lacked specific data on indicators of breeding (i.e., brood patch or enlarged cloaca) and detailed collection dates.
It is possible that our OK samples were migrating birds, but more sampling is required to ascertain whether white-crowned sparrows breed in OK and have this interesting partition of genetic structure by ecosystem, especially considering the possibility of notable range shifts as have been observed in Z. l. pugetensis in the last 35 years (Hunn & Beaudette, 2014).
Ecosite type also corresponded to morphological differences in white-crowned sparrow, albeit not by subspecies phenotype specifically ( Figure 3). Birds found in riparian deciduous forests were larger than alpine coniferous birds (Figure 4). The same phenomenon of larger body size in deciduous compared with coniferous ecosystems was observed in male pied flycatchers (Ficedula hypoleuca) resulting from competition (Lundberg et al., 1981). The clustering of white-crowned sparrows both genetically and morphologically into coniferous or deciduous forests is surprising. Most descriptions of white-crowned sparrow habitat include thickets and streamside shrubs (Z. l. gambelii), meadows of high sagebrush near the conifer timberline (Z. l. oriantha), and early seral coniferous forest (Z. l. pugetensis) (Dunn et al., 1995;Morton, 2002;Rivers et al., 2019) rather than along deciduous forest edge as we observed for some of our samples. Since aspen understory and canopy support greater invertebrate abundance than conifers (Rumble et al., 2001), it is possible that larger body size of white-crowned sparrows in deciduous forests is due to a greater abundance of high-quality food. Similar genetic patterns were found in black-capped chickadees, Poecile atricapillus, from the same area where patterns of genetic differentiation correspond to riparian poplar species and their hybrids (Adams & Burg, 2015). Body size may have also covaried with elevation, since 30% of our conifer forest samples were also from high elevations (>1,000 m). The greater energy costs of living in harsher climates of the high elevation sites could explain the smaller body size of sparrows sampled in those areas (Bonier et al., 2014). Although morphological patterns are congruent with genetic patterns, we cannot rule out the role of plasticity influencing morphological patterns. Therefore, we suggest that future studies should incorporate next-generation sequencing techniques to further explore the relationship between morphological variation and genetic variation.  (Table 1). Chickadees (Poecile spp.) from the same area show unusual genetic patterns. Hybrids of mountain (Poecile gambeli) and black-capped chickadees are more abundant in the FTSJ area (Grava et al., 2012), and populations of black-capped chickadees show high levels of population genetic structure (Adams & Burg, 2015).

| Species distribution models and landscape resistance
SDMs and LCCs showed that whereas elevation is an important contributor to modeling suitable habitat for three white-crowned sparrow subspecies, the Rocky Mountain range is a porous landscape barrier, especially for Z. l. gambelii. There is a lack of conspicuous barriers and a high degree of habitat suitability encompassing our Z. l. gambelii sample sites and much of the rest of North America ( Figure 5) which could explain why geographic distance (IBD) instead of landscape resistance (IBR) is the most prevalent barrier to dispersal (Table 4). IBD also explains most of the genetic distance between Z. l. oriantha populations, but unlike in Z. l. gambelii, IBR had small but significant importance at various spatial scales in Z. l. oriantha (Table 4). This difference may exist because Z. l. gambelii is a habitat generalist while Z. l. oriantha is more of a habitat specialist (Hahn et al., 1995;Maney et al., 1999;Wingfield et al., 1996).
Z. l. oriantha show robust responses to temperature and snowpack conditions in montane ecosystems, have strong site-fidelity, and require specific conditions for nesting (Wingfield et al., 2003).
Z. l. gambelii are "spatial opportunists" and show low site-fidelity (Hahn et al., 1995). This phenomenon is reflected in sympatric species of butterflies, in which the genetic structure of the generalist species was unaffected by the landscape matrix; however, the specialist species were highly sensitive to fine-scale ecosystem features (Engler et al., 2014).
In contrast to Z. l. oriantha and Z. l. pugetensis, the wider distri-  (Richardson et al., 2014). Our results highlight the importance of ecosystem in contemporary genetic differentiation (Figures 3 and 4), and LCC models show each subspecies' distribution is affected differently by environmental variables (Figure 5). Temperature is the strongest contributor for the Z. l. gambelii SDM, precipitation for Z. l. pugetensis, and elevation is important in the distribution models of all three subspecies, especially Z. l. oriantha. Further research is required to elucidate whether the current subspecies distributions are the result of historical range expansion (Banks, 1964;Rand, 1948), habitat fragmentation (Weir & Schluter, 2004), or recent divergence due to local adaptation of this ubiquitous species (Richardson et al., 2014).

| Genetic variation among subspecies
Genetic partitioning among subspecies of white-crowned sparrows has been previously established (MacDougall-Shackleton & MacDougall-Shackleton, 2001), and results indicate genetic partitioning among the subspecies investigated in this study provide further support for these patterns. Control region sequences indicate that Z. l. gambelii are distinct from Z. l. oriantha populations, although genetic variation between subspecies is relatively subtle. Only two control region haplotypes are shared between these two subspecies, and these haplotypes are found in southern Canada in the putative contact zone between these two subspecies. Z. l. pugetensis is also distinct from Z. l. gambelii, but shares several haplotypes with Z. l. oriantha. Microsatellite genetic patterns also indicate genetic differentiation between the three subspecies, and again individuals with admixed genotypes are found in southern Alberta where subspecies come into contact. Although other studies have found limited support for genetic distinctiveness among subspecies (Ball & Avise, 1992;Zink et al., 2000), genetic patterns do not always follow subspecies boundaries (Zink, 2004). The prevalence of genetic partitioning among the subspecies investigated in this study is further evidence of how complex genetic patterns are for this species and indicates that current taxonomy should continue to recognize the existing subspecies.

| CON CLUS ION
In conclusion, it is important to understand all the influences and barriers on population genetics in montane systems as they are an important ecological zone for biodiversity (Fjeldsâ et al., 2011). The western extent of the white-crowned sparrow's distribution is genetically structured not by subspecies, but by the combination of distance, mountains, and ecotype. White-crowned sparrows showed both genetic differentiation and morphological differentiation in coniferous versus deciduous ecosystems, a phenomenon of habitat partitioning increasingly observed in other population genetics studies (Grava et al., 2012;Jensen et al., 2019;Lundberg et al., 1981;Porlier et al., 2012). Using a widespread and ubiquitous avian species can contribute to an understanding of barriers that could be affecting the gene flow in species of conservation concern which are more difficult to access. For populations of threatened species isolated after habitat fragmentation and degradation, using a model species to understand the factors affecting dispersal and gene flow can be critical for conservation and management decisions (Dubey et al., 2011).

ACK N OWLED G M ENTS
We

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data are archived in Dryad https://doi.org/10.5061/dryad.zw3r2 2882.