Agriculture creates subtle genetic structure among migratory and nonmigratory populations of burrowing owls throughout North America

Abstract Population structure across a species distribution primarily reflects historical, ecological, and evolutionary processes. However, large‐scale contemporaneous changes in land use have the potential to create changes in habitat quality and thereby cause changes in gene flow, population structure, and distributions. As such, land‐use changes in one portion of a species range may explain declines in other portions of their range. For example, many burrowing owl populations have declined or become extirpated near the northern edge of the species' breeding distribution during the second half of the 20th century. In the same period, large extensions of thornscrub were converted to irrigated agriculture in northwestern Mexico. These irrigated areas may now support the highest densities of burrowing owls in North America. We tested the hypothesis that burrowing owls that colonized this recently created owl habitat in northwestern Mexico originated from declining migratory populations from the northern portion of the species' range (migration‐driven breeding dispersal whereby long‐distance migrants from Canada and the United States became year‐round residents in the newly created irrigated agriculture areas in Mexico). We used 10 novel microsatellite markers to genotype 1,560 owls from 36 study locations in Canada, Mexico, and the United States. We found that burrowing owl populations are practically panmictic throughout the entire North American breeding range. However, an analysis of molecular variance provided some evidence that burrowing owl populations in northwestern Mexico and Canada together are more genetically differentiated from the rest of the populations in the breeding range, lending some support to our migration‐driven breeding dispersal hypothesis. We found evidence of subtle genetic differentiation associated with irrigated agricultural areas in southern Sonora and Sinaloa in northwestern Mexico. Our results suggest that land use can produce location‐specific population dynamics leading to subtle genetic structure even in the absence of dispersal barriers.


| INTRODUC TI ON
Understanding ecological and evolutionary dynamics of a species at the edges of its distribution can help unveil the mechanisms that limit abundance throughout a species' entire geographic range (Holt & Keitt, 2005). In this regard, ecological theory and empirical evidence support the idea that species tend to be less abundant and more prone to local population extinction at the periphery of their geographic ranges (Gaston, 2003). Populations at the edge of a species' distribution may be maintained by dispersal and recolonization from interior populations (Curnutt, Pimm, & Maurer, 1996).
This scenario whereby populations on the periphery are repeatedly "rescued" (Brown & Kodric-Brown, 1977) by interior populations may be particularly important for species of conservation concern.
Understanding the processes by which peripheral populations are maintained in those species is important for designing effective recovery efforts. For example, populations of the western burrowing owl (Athene cunicularia hypugaea) have been extirpated from some areas and are rare and declining in other areas near the northern edge of their breeding distribution (Clayton & Schmutz, 1999;Macías-Duarte & Conway, 2015;Skeel, Keith, & Palaschuk, 2001;Wellicome & Holroyd, 2001). The breeding range of the western burrowing owl (burrowing owl hereafter) historically comprised semiarid grasslands from southern Canada to central Mexico (Poulin, Todd, Haug, Millsap, & Martell, 2020). Hypotheses to explain population declines in the northern portion of their range include local mechanisms such as conversion of grassland to dryland farming in the northern Great Plains, extirpation of black-tailed prairie dogs, toxicological effects of pesticides, collisions with vehicles, and annual dispersal (Clayton & Schmutz, 1999;Desmond, Savidge, & Eskridge, 2000;Duxbury, 2004;Klute et al., 2003;Poulin et al., 2020). All these hypotheses seem insufficient to explain the extent of burrowing owl population declines observed in the northern portion of their breeding range because many areas with seemingly suitable habitat remain unoccupied. Nevertheless, the highest breeding densities occur in the southern portion of the burrowing owl's breeding range (in Imperial Valley, California; DeSante, Ruhlen, Rosenberg & Haley, 2004;Sauer et al., 2017). In addition, densities of breeding burrowing owls in the coastal plains of Sonora and Sinaloa may be just as high as those in southeastern California (Macías-Duarte, 2011). These high densities of burrowing owls in the southern portions of the species' range are all in arid desert areas that have been converted to irrigated agriculture. High densities of breeding burrowing owls in this portion of their range is a recent phenomenon; more than 1.5 million hectares of coastal thornscrub and tropical dry forest in Sonora and Sinaloa were converted to irrigated farmland in the last 60 years (Anonymous, 1994;Rohwer, Grason, & Navarro-Sigüenza, 2015). This redistribution of burrowing owls (the breeding range contracting in the north and expanding in the south) poses interesting questions about the mechanisms that shape and maintain the geographic range of the species especially given that many other birds in North America are showing opposite trends (ranges shifting northward) in response to climate change (La Sorte & Thompson, 2007). In this paper, we propose and test the hypothesis that the contraction at the northern periphery of the burrowing owl's range and the expansion in the southern portion of their range may be directly related.
Most breeding populations of burrowing owls in North America exhibit partial migration (where some individuals migrate and some do not), but northern populations in the Great Plains are 100% migratory (Poulin et al., 2020). James (1992) speculated that burrowing owls have a leap-frog migration pattern (negative correlation between breeding latitude and wintering latitude across populations).
In addition, most burrowing owls that breed in the northern portion of the breeding range appear to spend their winters in southern Texas and central Mexico (Duxbury, 2004;Holroyd, Trefry, & Duxbury, 2010). We tested the hypothesis that burrowing owls that once migrated annually from northern portions of their breeding range to central Mexico became resident breeders in these newly created irrigated agricultural areas, contributing to both population declines in the north and population increases in the south. Birds breeding within what was formerly their wintering grounds (migrants becoming year-round residents) has been contemporarily observed in at least 3 other species (Sutherland, 1998). However, numerous phylogenetic analyses infer that these migratory dropoffs have been common through the evolutionary histories of migratory birds and are drivers of diversification and speciation (e.g., Gómez-Bahamón et al., 2020;Rolland, Jiguet, Jønsson, Condamine, & Morlon, 2014;Voelker & Light, 2011;Winger, Barker, & Ree, 2014).
Testing this hypothesis requires inferring patterns of breeding dispersal (movement between two breeding attempts; Greenwood & Harvey, 1982) among populations throughout the burrowing owl breeding range. We used genetic markers to infer the patterns of gene flow produced by breeding dispersal by measuring genetic differentiation among migratory and nonmigratory populations throughout North America. We tested 3 predictions of our hypothesis that infer patterns of genetic variation produced by gene flow from northern migratory (declining) populations to southern populations within irrigated agricultural areas. First, our hypothesis predicts that genetic differentiation between a northern migratory population and a southern agricultural population will be lower than the expected genetic differentiation predicted by the geographic distance between the 2 populations. This prediction assumes an isolation-by-distance pattern (Wright, 1943), where populations further apart geographically are more genetically differentiated (however, subtle) than populations closer to each other due to differences in Athene cunicularia hypugaea , dispersal, DNA microsatellites, gene flow, genetic differentiation, irrigated agriculture, migration, population genetics frequency of dispersal. Second, our hypothesis predicts that northern migratory populations and southern agricultural populations together are genetically similar enough to be differentiated from the rest of the breeding populations within the burrowing owl breeding range. This prediction can be tested via a significance test of the two-group classification of burrowing owl populations mentioned above to explain overall genetic variation. We can use an assignment test to test a third prediction. Assignment tests use individual genotypes to estimate the probability of membership of each individual genotype to predefined clusters of individuals. In this regard, our hypothesis predicts that southern agricultural populations will have more individual owls with probabilities of membership similar to those found in individuals from northern migratory populations (in areas where owls are declining) compared to the nonagricultural populations in the southern part of the species' range. We used DNA samples from owls throughout their North American breeding range to test these 3 predictions.

| Study area
We obtained DNA samples from 1,560 breeding burrowing owls from 36 locations ('study locations' hereafter) in Canada, Mexico, and the United States ( Figure 1, Table 1). To test our predictions,  Table 1. Yellow labels denote northern declining migratory populations and blue labels denote southern agricultural populations. The gray area denotes the breeding distribution of the burrowing owl(after Poulin et al., 2020) owl, but systematic regional declines have been most evident in Alberta, Saskatchewan, North Dakota, and South Dakota, where the species is close to extirpation (owls have been extirpated from Manitoba and British Columbia). Therefore, we only defined Alberta (ALB), Saskatchewan (SAK), and Grand River-Little Missouri National Grasslands (GRL) as northern study locations with declining migratory breeding populations ("northern study locations" hereafter) ( Table 1).

| Sample collection
We trapped burrowing owls during the summers of 2004-2009. We trapped burrowing owls using push-door tramps (Winchell, 1999) set at the entrance of nest burrows and bownet traps (Bub, 1991) set near nest burrows. None of the 1,560 birds that we included in our analysis were closely related (i.e., a parent and its offspring, or >1 juvenile from the same nest burrow). Our primary source of genomic DNA was blood. We obtained ~50 μl of blood through a venipuncture of the brachial vein. We also used flight and/or body feathers occasionally as a source of genomic DNA when we could not withdraw a blood sample.

| Genotyping
We used 10 microsatellite markers developed specifically for this study (

| Data analysis
We used MS Excel© macro GENALEX 3.6 (Peakall & Smouse, 2006) to calculate standard descriptive statistics of genetic diversity of burrowing owls at each of our study locations, including observed heterozygosity, expected heterozygosity, and fixation index F. We also used program ARLEQUIN 3.1.1 (Excoffier, 2006) to estimate the Weir and Cockerham's F ST (θ, Weir & Cockerham, 1984) for all populations.
We computed actual differentiation D (measure of differentiation between populations independent of gene diversity) (Jost, 2008) to test our prediction that gene flow between declining migratory populations in the north and populations in southern agricultural areas would disrupt an otherwise apparent isolation-by-distance relationship. We used software SMOGD (Crawford, 2010)to compute actual differentiation D. We used D as our measure of population-pairwise genetic differentiation because F ST does not adequately measure genetic differentiation when within-population allelic diversity is high (Jost, 2008). D ranges from 0 to 1, corresponding to complete similarity to complete differentiation. We performed a Mantel test (Mantel, 1967) to test our assumption of the existence of an isolation-by-distance pattern (i.e., that the genetic differentiation between 2 populations is positively correlated to the geographic distance that separates those populations  (Michalakis & Excoffier, 1996).
The former measure assumes the stepwise mutation model (Ohta & Kimura, 1973), which is appropriate for microsatellite loci. We We conducted an assignment test as implemented by the program STRUCTURE (Hubisz, Falush, Stephens, & Pritchard, 2009;Pritchard, Stephens, & Donnelly, 2000) to test our prediction that southern agricultural study locations will have more individual owls with probabilities of membership similar to those found in individuals from declining populations in the north compared to the nonagricultural study locations in the southern part of the species range. STRUCTURE 2.3.3 implements an algorithm suited to infer weak population structure (Hubisz et al., 2009). We used study locations as prior information to assist the inference of population structure (Hubisz et al., 2009). We performed 10 runs for each K = 1, 2, … 10. Each run consisted of a burn-in period of 50,000 Markov Chain Monte Carlo repetitions followed by 50,000 repetitions to sample from the posterior distribution of K. We estimated L(K) for each K from correlated allele frequencies and an admixture model. This approach is superior at detecting subtle genetic structure when population differentiation is low compared to the use of uncorrelated allele frequencies and a nonadmixture model (Falush, Stephens, & Pritchard, 2003 it performs better in detecting population genetic structure than L(K) (Evanno, Regnaut, & Goudet, 2005). Therefore, actual number of populations is revealed by the value of K with the highest value of ΔK. We used program CLUMPP (Jakobsson & Rosenberg, 2007) to calculate the posterior probabilities of membership of each individual owl to each of the K populations from our multiple runs in STRUCTURE. Population acronyms are shown in Table 1. Study location acronyms with (*) and ( †) denote southern agricultural populations and northern declining migratory populations, respectively. Different shades denote population membership based on >50% of the posterior probability of membership from the program STRUCTURE (see Figure 3).

| RE SULTS
Burrowing owls exhibited high levels of genetic diversity (Table 2) with relatively low variation among study locations. Per-locus av-  additional AMOVAs based on F ST was significant for both the standard AMOVA and the weight-averaged AMOVA over all loci (Table 3), which is precisely the AMOVA that included the nearest 3 study locations (CHI, JAN, and TUC) within Group 1.
STRUCTURE revealed a genetic structure consisting of 3 populations in the burrowing owl in spite of the low levels of genetic differentiation among study locations shown by F ST and D statistics. Mean log-likelihood of the observed genotypic data and ΔK was highest at K = 3 (indicating 3 distinct populations, see Figure S1). The posterior probabilities of membership of each of our 1,560 individual owls assigned to these putative populations (see Figure S2)

| D ISCUSS I ON
Burrowing owl populations in North America have low levels of differentiation as shown by F ST and D statistics and, in that regard, our results corroborate a previous study that also reported low levels of genetic differentiation for the western burrowing owl (Korfanta, McDonald, & Glenn, 2005 ). Korfanta et al. (2005) Korfanta et al., 2005). Therefore, burrowing owls clearly have low genetic differentiation among populations that extends throughout the entire breeding range in North America (including populations in Canada and Mexico), which is not surprising for a migratory bird.
However, a major assumption for our 3 predictions is that burrowing owl populations had at least some subtle genetic structure before the development of the agricultural valleys in southwestern United States and northwestern Mexico. This low genetic population differentiation throughout the burrowing owl breeding range hindered our ability to rigorously test the migration-driven breeding dispersal hypothesis. Despite the challenges associated with the minimal genetic structure, we did detect some tentative support for the migration-driven breeding dispersal hypothesis.
The 10 DNA microsatellite loci we used may have been too few to detect subtle genetic structure using the F ST index. Hence, the use of our genetic markers to detect past and current patterns of breeding dispersal is imperfect. We used numerous analytical methods and algorithms that made use of the individuals' genotypic data  did not provide support of the hypothesis. However, measures of allele size have been criticized for having large sampling errors and low efficiency in reconstructing simulated phylogenies (Takezaki & Nei, 1996). In addition, the lack of statistical significance in 6 of the 7 additional AMOVAs (Table 3) Hull et al., 2010). Low levels of genetic differentiation among populations of burrowing owls are highly relevant for burrowing owl conservation and restoration programs everywhere in North America.
Low genetic differentiation among our 36 study locations from Canada to central Mexico provides further evidence that burrowing owls are a large panmictic population across the species' breeding range. Reintroduction programs may be able to use individuals from populations throughout western North America without substantially compromising genetic variation for local adaptation. Low genetic differentiation, presumably caused by continent-wide breeding dispersal, also means that population trends in a given location may be caused by changes in demographic processes (e.g., fecundity, mortality, and emigration) in other portions of the species' range.
Therefore, population declines in the northern edge of the species' breeding distribution may reflect either declines in immigration from more interior populations, low local recruitment, or both.

ACK N OWLED G M ENTS
We complied with Canadian, Mexican, and U. S. regulations for bird handling, blood and feather collection, and import-export across international boundaries. We thank the U.S. Department of Defense

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