Demographic history shapes North American gray wolf genomic diversity and informs species' conservation

Effective population size estimates are critical information needed for evolutionary predictions and conservation decisions. This is particularly true for species with social factors that restrict access to breeding or experience repeated fluctuations in population size across generations. We investigated the genomic estimates of effective population size along with diversity, subdivision, and inbreeding from 162,109 minimally filtered and 81,595 statistically neutral and unlinked SNPs genotyped in 437 grey wolf samples from North America collected between 1986 and 2021. We found genetic structure across North America, represented by three distinct demographic histories of western, central, and eastern regions of the continent. Further, grey wolves in the northern Rocky Mountains have lower genomic diversity than wolves of the western Great Lakes and have declined over time. Effective population size estimates revealed the historical signatures of continental efforts of predator extermination, despite a quarter century of recovery efforts. We are the first to provide molecular estimates of effective population size across distinct grey wolf populations in North America, which ranged between Ne ~ 275 and 3050 since early 1980s. We provide data that inform managers regarding the status and importance of effective population size estimates for grey wolf conservation, which are on average 5.2–9.3% of census estimates for this species. We show that while grey wolves fall above minimum effective population sizes needed to avoid extinction due to inbreeding depression in the short term, they are below sizes predicted to be necessary to avoid long‐term risk of extinction.


| INTRODUC TI ON
The theory of the effective population size (N e ) was originally developed by Sewall Wright (1943Wright ( , 1965) ) to provide a means for comparing structure across seemingly disparate populations to result in an estimate that represents an idealized population of randomly mating individuals (Crow & Kimura, 1970).Thus, social organization and non-random breeding will impact the distribution of genotypes over geographic space and concomitantly N e estimates.Any factor that results in deviations from random breeding (e.g.social factors, breeding strategies, site availability) or changes in population size across generations will result in an effective population size estimate that is a fraction of the census size (N) (Charlesworth & Willis, 2009;Clutton-Brock, 2016;Hedrick & Kalinowski, 2000;Keller & Reeve, 1994).For species with high reproductive skew and social structures that repress reproduction in subdominant ranks, the effective population size estimate inferred from sex ratios, dispersal or migration rates, number of reproductive individuals, or genetic assessments is critical information needed for evolutionary predictions (Lanfear et al., 2014;Wang et al., 2016).
Population sizes fluctuate over time, either through natural process or due to anthropogenic activity such as wildlife management (Rowe & Beebee, 2004).Any reduction in size, compounded with isolation, will erode genetic variation via random genetic drift to a degree that depends on the severity and duration of these bottlenecks (Fisher, 1958).Without inter-population connectivity, the only process that naturally introduces new variation into the gene pool is de novo mutations.New mutations are more likely to quickly drift to fixation in isolated small populations, resulting in continuing low levels of genetic diversity (Coyne et al., 1997;Fisher, 1930;Wade & Goodnight, 1998;Wright, 1931).The potential for a population to respond to evolutionary challenges deteriorates as genomic variation dwindles, thereby limiting adaptive outcomes (Allendorf, 2016;Frankham, 2005;Hoffmann et al., 2017;López-Cortegano et al., 2019).Anthropogenic effects that reduce population size and impact life history events central to individual-level fitness (e.g.reproduction, dispersal) are well known to degrade genomic variation and adaptive potential (Allendorf et al., 2008;Coltman, 2008;Frankel & Soulé, 1981;Frankham, 2005;Reed & Frankham, 2003).
In their recent evolutionary history, grey wolves (Canis lupus) in North America have been eradicated from much of their southern continental range through federal and state programmes first implemented during the mid-19th century.These programmes were highly effective and by the late 1950s had exterminated the wolf from the conterminous United States except for a few individuals on Isle Royale National Park in Lake Superior (Minnesota) and a few hundred individuals in northeastern mainland Minnesota (Boitani, 2003;Franzmann & Schwartz, 1997;Kolenosky & Standfield, 1975;Parker, 1995;Peterson, 1955;U.S. Fish and Wildlife Service, 1992;Young & Goldman, 1944).In the face of a near total elimination, coupled with social structure of the species and removal of dispersers, there was a growing concern regarding the future survival of the grey wolf species which led to the translocation of grey wolves to Yellowstone National Park (YNP) and central Idaho (Adams et al., 2008;Brainerd et al., 2008;Rick et al., 2017;Treves et al., 2016).A targeted study of wolves living within YNP reported a significantly smaller effective population size than the censused population (vonHoldt et al., 2008), emphasizing the critical role of population connectivity to combat genetic drift, inbreeding, and erosion of heterozygosity (Allendorf et al., 2008;Gese & Mech, 1991;Jedrzejewski et al., 2005;Mech & Boitani, 2003;von-Holdt et al., 2008).
In the United States, grey wolves are managed as three populations with distinct demographic histories: northern Rocky Mountains, the western Great Lakes, and southwestern (explicitly  (Refsnider, 2009).The Timber Wolf Recovery Plan further considered the historic range to Minnesota eastward to Maine and south to the northern portion of Florida (Refsnider, 2009;U.S. Fish and Wildlife Service, 1992;Wisconsin, 1989;Wydeven et al., 2009).The southwestern population that encompasses the endangered Mexican grey wolf subspecies was not included in this study.
Effective in January 2021, the U.S. Fish and Wildlife Service (FWS) delisted grey wolves (excluding the Mexican wolf subspecies) everywhere in the lower 48 United States (final rule 85 FR 69778).
By February 2022, ESA protections were restored for all grey wolves in the lower 48 United States except for the wolves of the northern Rocky Mountain region, where they remain under state-level management.The delisting decision relied in part on the lack of information from FWS that the western Great Lakes population could indeed be self-sustaining without federal protection.By January 2023, the Circuit Mediator issued an order for a scientific review of grey wolf status review to be conducted.
Our goal was to assess the temporal and spatial variations in genetic signatures over the recent decades of grey wolf protections and recovery across portions of North America and provide information to consider for long-term viability of grey wolves as it pertains to their ESA listing status in the United States.We conducted this genomic surveillance across the North American continent to showcase how demography and genomic signatures are intertwined.
This assessment provides a contemporary assessment of genetic parameters important to genomic viability across geographic and regulatory scales for integration into conservation goals for a social carnivore species.

| Sample collection and genomic library construction
We obtained archived blood or tissue samples collected from 482 grey wolves across their continental range in North America (Canada = 91, USA = 391) from state and federal partners, local trappers, and private genetic collections (Figure 1a; Table S1).
Locations of sample origins varied, from regional identification to counties, parks, or states and provinces.We partitioned samples into two levels of geographic resolution, regional and U.S.-managed populations.For the U.S.-managed populations, we define the 'northern Rocky Mountains' (abbreviated as RM) as composed of samples that originated from California, Idaho, Montana, Washington, and Wyoming.We define Michigan, Minnesota, and Wisconsin to compose the 'western Great Lakes' (abbreviated as GL).
We extracted genomic DNA following manufacturer's protocol (Qiagen DNeasy Blood and Tissue kit).We used the Qubit fluorometer system for DNA quantification to standardize the input amount for use in the modified restriction-site-associated DNA sequencing (RADseq) capture protocol (Ali et al., 2015).Briefly, we digested genomic DNA with SbfI with a subsequent ligation of unique 8-bp barcoded biotinylated adapters to permit the pooling of 48 DNA samples into a single library.We randomly sheared each library to 400 bp in a Covaris LE220 followed by an enrichment for the adapter-ligated fragments using a Dynabeads M-280 streptavidin binding assay.We then prepared the enriched libraries for pairedend (2 × 150 nt) Illumina NovaSeq 6000 sequencing at Princeton University's Lewis-Sigler Genomics Institute core facility using the NEBnext Ultra II DNA Library Prep Kit (New England Biolabs).For any step of purifying or size selection of DNA, we used Agencourt AMPure XP magnetic beads (Beckman Coulter).

| Bioinformatic processing
We retained sequence read pairs that contained both our known unique barcodes and remnant SbfI recognition site, which were processed in STACKS v2.6 (Catchen et al., 2013;Rochette et al., 2019).We used the process_radtags module to rescue our barcoded reads with a 2 bp mismatch and excluded reads with a quality score < 10.We next removed PCR duplicates in the clone_filter module followed by mapping to the reference dog genome CanFam3.1 assembly (Lindblad-Toh et al., 2005) using bwa-mem (Li, 2013).We also included the Y chromosome (KP081776.1;Li et al., 2013) with the CanFam3.1 reference assembly.After alignment, we excluded mapped reads with MAPQ <20 and then converted the SAM files to BAM format in Samtools v0.1.18(Li et al., 2009).We implemented the gstacks and populations modules in STACKS v2 with an increase in the minimum significance threshold in gstacks and used the maximumlikelihood marukilow model that incorporates uncertainties for lowcoverage data (-vt-alpha and -gt-alpha with p = .01).We additionally used the flag -r 60 to retain only newly annotated sites found in at least 60% of the samples in the catalogue.In VCFtools v0.1.17(Danecek et al., 2011), we estimated the pre-filtered sequence coverage and then subsequently filtered loci to exclude singleton and private doubleton alleles, removed loci with more than 90% missing data across all samples, and excluded individuals with more than 30% missing data.We removed loci with a minor allele frequency (MAF < 0.03) and required at least an 80% genotyping rate per locus (-geno 0.2) in PLINK v1.90b3i (Chang et al., 2015).
We used VCFtools for individual-level metrics of heterozygosity (observed, H O ; expected, H E ) and the two-sample Kolmogorov-Smirnov to test for statistical differences in data distributions and correlations in R (R Core Team, 2022).We then utilized the populations module in STACKS v2 to identify alleles private to each canid lineage.We further conducted a rarefaction method for private allele richness per locus while controlling for sample size variation in the number of genomes sampled in the programme ADZE (Szpiech et al., 2008) with the parameter G of sample size set to 100.

| Sex inference from sequence coverage of the Y chromosome
As we included the Y chromosome (KP081776.1;Li et al., 2013) with the CanFam3.1 reference assembly for read alignment, we used ttests and the two-sample Kolmogorov-Smirnov to determine the sequence coverage differences between the sexes.This provided us an opportunity to establish a threshold of Y-specific sequence coverage to infer sex, with females inferred from falling below the threshold and males above.We then repeated analyses independently for each sex to explore the impact of sex-biased demography.

| Population structure and differentiation
For demographic analyses, we constructed a statistically neutral and unlinked dataset of SNPs by excluding sites within 50-SNP windows that exceeded genotype correlations of r = .2(-indep-pairwise 50 5 0.2; a proxy for linkage disequilibrium or LD) and SNPs that significantly deviated from Hardy-Weinberg Equilibrium (HWE) with the argument -hwe 0.001.We conducted both non-model and model-based clustering analyses.We completed the former as a principal component analysis (PCA) in FlashPCA v2. 1 (Abraham et al., 2017) and the latter with an unsupervised maximum-likelihood framework with Admixture (Alexander et al., 2009).We analysed the fit of two to 10 partitions (K) with the cross-validation error (cv) flag.
We also estimated inter-group pairwise genetic differentiation as Weir and Cockerham's F ST in VCFtools v0.1.17.We reported average F ST across the genome (autosomes and X chromosome combined).

| Inbreeding estimates from autozygosity
We analysed the minimally filtered SNP set separately for loci on the autosomes and X chromosome.These loci represented a total F I G U R E 1 Population genetic structure of 437 grey wolves from (a) North American populations genotyped at 81,595 statistically neutral and unlinked SNPs inferred from (b) principal component analysis (axes rotated to show geographic correspondence); and (c) a maximumlikelihood approach for three and nine partitions (map credit: Free Vector Maps WRLD-NA-01-0007).(d) Rarefaction of allelic richness and private alleles for each major geographic region of grey wolves (see Table S1).(Table S1).To detect autozygosity from runs of homozygosity (ROH), we used the following parameters for low-coverage data: homozygdensity 50, homozyg-gap 1000, homozyg-kb 300, homozyg-snp 50, homozyg-window-het 4, homozyg-window-missing 5, homozygwindow-snp 50, and homozyg-window-threshold 0.05 (Ceballos et al., 2018).We converted the ROH segments to an individual-level inbreeding coefficient (F ROH ) following Taboada et al. (2014): where L ROH is the length of an ROH segment in an individual.

| Effective population size estimates
We estimated effective population (N e ) sizes and focused on recent (past 200 generations) estimations as presumed to be more accurate.
Effective population size estimates extrapolate population parameters from genetic diversity metrics.Although dispersal and translocation events are known, the collection of genetic variation is the core of such inference and is bounded by how a population is defined in time and space.Here, we implemented the algorithm in GONE (Santiago et al., 2020), which is an LD-based method that accounts for drift (i.e.finite census size) and makes use of recombination rates but is influenced by both population structure and admixture.GONE leverages a genetic algorithm from Mitchell (1998) to search across sequences of possible historical effective population sizes that best explain the spectrum of observed LD values to minimize the sum of squares of the differences between observed and expected allelic covariances.We assumed unphased data, no MAF pruning, a maximum of 50,000 SNPs considered per chromosome, and ignored pairs of SNPs with recombination rate over 0.05, as recommended for the software.A constant rate of recombination of 1 cM per Mb was assumed across the genome.We estimated N e sizes at two levels: each major geographic region and population designations for management implications in the United States.However, resulting estimates for the wolf populations in Canada should be interpreted with caution given our limited genotype surveillance across the region.We estimated N e from autosomal SNP data and translated generations into years using 4 years per generation as the unit of time (Mech et al., 2016;vonHoldt et al., 2008).We believed that only the minimally filtered RADseq data (i.e.missingness and MAF) was appropriate for these estimates (Beichman et al., 2017).
Finally, we were conservative when interpreting 'present-day effective population size' as the most recent four generations for N e are considered a single analytical block by GONE.Hence, we used the N e average of generations 1-8 to avoid biases from any lingering artefact in generations 1-4 (Novo et al., 2023).We also focus on reporting the results of the last 50 generations (approximately 200 years) as that is most pertinent to the recent population demography and conservation considerations.
We then assessed how well the effective population size estimates explain the expected decay in heterozygosity using the formula when t = 8:

| Admixture is part of the history of the western Great Lakes grey wolf population
We rediscovered SNPs with the addition of BAM files from previ-  S1b).The grey wolves in the Great Lakes region are known to have a history of admixture with both coyotes and eastern wolves (Heppenheimer et al., 2018;vonHoldt et al., 2011).The predominant signal described to date is that Great Lakes region grey wolves have partial coyote ancestry with grey wolves of southeastern Ontario carrying more partial ancestries of eastern wolves.These were merged with the BAM files from the population of northern Rocky Mountains and western Great Lakes samples to explore the impact of coyote and eastern wolf admixture on grey wolf genetic estimates.We followed the same analysis and filtering methods as described above to obtain a statistically unlinked and neutral set of SNP loci.We conducted an unsupervised assignment analysis for K = 2-10 in ADMIXTURE and complemented with genetic differentiation (F ST ) estimates using VCFtools v0.1.17.

| Reliable inferences from reduced representation low-coverage population-level genotype data
Population genomic studies can leverage the affordable technologies of reduced representation data collection methods, such as RADseq, to collect genotype data from hundreds or thousands of individuals.The drawbacks are obvious in terms of missing rare alleles or allele dropout rates due to the nature of the library preparation.
Thus, studies have assessed the biases and challenges of lowcoverage data (3-6×) compared to whole-genome sequence (WGS) and found that the former can be equally informative with careful adjustments to methods and inferences (Ceballos et al., 2018;Duntsch et al., 2021).It is known that some population metrics like ROH are expected to be biased.For example, low-coverage data likely underestimate the frequency of small and overestimate larger ROH fragments (Lavanchy & Goudet, 2023).

| RE SULTS
We sequenced 482 grey wolf samples from North America, collected between 1986 and 2021 when known, with an average fold sequence coverage of 7.3 (±3.4) to discover 1,099,764 raw, RAD loci that passed our STACKS filtering parameters but prior to populationlevel filtering (Table S1).We excluded 45 wolves due to high (>20%) When we presumed these samples having correct sex inference, the average sequence depth on the Y chromosome was significantly enriched in males (females = 3406.9,males = 25587.Kolmogorov-Smirnov D = 1.0, p < 10 −16 ) (Figure S1b).

| Grey wolves are genetically and geographically structured across North America
We presented two levels of genetic structure across the North American continent that reflect the geographic assignment probabilities for two cluster analyses: the PCA (K = 3) and the best supported partition from maximum-likelihood inference (K = 9) (Figure 1b,c;  S1).The other two clusters represent northern Quebec and the shared demography of Ontario and the western Great Lakes population (Table S1).Out of these four geographic groupings, we found that only two groups carried private alleles (western USA and southwestern Canada, n = 332; Great Lakes and Ontario, n = 6801) out of 162,109 SNPs.A rarefaction analysis mirrors the demographic history of each, with the Great Lakes and Ontario regional group showing the highest level of allele richness and mean number of private alleles per locus controlled for sample size differences (Figure 1d), likely due to their known history of coyote and eastern wolf admixture (Koblmüller et al., 2009;vonHoldt et al., 2016).had the lowest levels of genetic differentiation while the highest was found between opposite coasts of the continent (western USA and southwestern Canada vs. northern Quebec, F ST = 0.084) (Table 1, Figure S3).We find that all genetic differentiation distributions are significantly distinct (Table S2).We assessed this metric for females and males separately for two geographic regions (western USA and southwestern Canada; Great Lakes and Ontario).While northern Rocky Mountain grey wolves showed variable levels of differentiation within the region (F ST genome = 0.0-0.13,X = 0.0-0.09),females were significantly higher levels of genome-wide differentiation to other females (female-female F ST = 0.052) than males (male-male F ST = 0.032, 1-tailed t-test of unequal variance p = .01207)(Figure S4a).In contrast, western Great Lakes grey wolves had much lower intra-region genetic differentiation (F ST genome = 0.0-0.03,X = 0.0-0.04),with no significant differences between males and females (F ST female-female = 0.017, male-male = 0.019, p = .3242)(Figure S4b).

| Genomic diversity and inbreeding coefficients are variable across continental North America
Northern Quebec grey wolves had the highest levels of observed and expected heterozygosity estimates (H O = 0.284), followed by  with the latter found to have significantly higher observed heterozygosity than expected (Table 2a).We further report the expected positive correlation between the number of autosomal ROH segments and inbreeding estimates (R = .77),with a weaker yet similar trend for the X chromosome (R = .44).Autosomal inbreeding levels were highest in the wolves of western USA and southwestern Canada (F ROH = 0.296), which were not significantly different from northern Canada (F ROH = 0.278) or northern Quebec (F ROH = 0.267).Wolves of the Great Lakes/Ontario (F ROH = 0.199) had significantly lower inbreeding levels (F ROH = 0.278) than the other geographic regions.

| The northern Rocky Mountain population is genetically distinct
To provide information relevant to ongoing management considerations and decisions, we partitioned the samples to analyse only those belonging to the populations identified in the United States,   H E : R = −.46,p = 1.2 × 10 −8 ) (Figure 2a).Although the WGL population shows a similar albeit weaker pattern of decline, there was no statistical significance (H O : R = −.08,p = .47;H E : R = −.12,p = .26)(Figure 2b).
Females in the northern Rocky Mountains population were significantly more differentiated from each other than males across the genome (mean F ST = 0.052 and 0.032, respectively; 1-tailed t-test of unequal variance p = .01207)and the X chromosome (F ST = 0.051 and 0.029; p = .0051)(Figure S4).This pattern was not found in the females of the western Great Lakes population (genome: F ST = 0.017 and 0.019; p = .3242;X chromosome: F ST = 0.016 and 0.012; p = .1876).
The northern Rocky Mountain grey wolves had significantly higher autosomal inbreeding coefficients compared to the western Great Lakes, which differences across the X chromosome were not significant (F ROH , autosomes: RM = 0.299, GL = 0.211, t = 8.5, df = 309.6,p = 8.67 × 10 −16 ; X chromosome: RM = 0.076, GL = 0.070, t = 0.8, df = 260.3,p = .4473)(Figure S5).The outlier inbreeding coefficients for western Great Lakes can be attributed to the small and isolated grey wolf population living in Isle Royale National Park.

| Population effective size estimates show the continental history of extermination and recovery
We inferred population effective sizes for the past 50 genera- the four regional genetic clusters that carried genetic distinction.
We estimated N e ranged between 63.0 and 3848.5 over the past 50 generations at a regional scale (Figure 3a; Table S3).Northern Canada had the highest historical size estimated at 3848.We found a significant positive relationship between regional effective population size and number of generations before present (Pearson's product-moment correlation R = .39,t = 7.3, df = 298, p = 2.03 × 10 −12 ) (Figure S6).When we restricted our analyses to the two populations, we found that the northern Rocky Mountains wolves, respectively (Table S3).
We further compared population estimates for the northern Rocky Mountains and western Great Lakes populations obtained from management, agency, and public reports between 1982 and 2015 (Table S4).(Figure 2a).When we use the estimated average effective population size N e = 141.7 for the northern Rocky Mountains during that time (Table S3), we estimate that the observed heterozygosity should decay by 0.032 to H O = 0.203, which is within the 95% confidence interval (Figure 2a).We found the same trend for the western Great Lakes (H O ~ 0.213 and 0.213 in 1988 and 2020, respectively), estimated to decay by 0.016 to H O = 0.197.

| Admixture with coyotes and eastern wolves is unique to the Great Lakes grey wolves
We created a second dataset that included western coyotes and eastern wolves to explore signatures of admixture in the grey wolves of the Great Lakes region.We discovered 163,314 genomic loci Mountains population, 184 from the western Great Lakes population, 74 western coyotes, and 28 eastern wolves).We also constructed a statistically neutral and unlinked dataset of 80,655 SNPs.
At the highest level of partition (K = 10), we found that grey wolves of the western Great Lakes population had the highest average (±sd) probability assignment to clusters of other Great Lakes grey wolves (Q = 0.64 ± 0.4) and <10% to any other wolf group (3.4 ± 0.1% assignments to eastern wolves; <2% to Rocky Mountain grey wolves), with minimal assignments to western coyotes (Q = 0.01 ± 0.1) (Figure S7; Table S5).Rocky Mountain grey wolves similarly formed their own cluster (Q > 0.97) with low, albeit detectable, partial assignments of Wyoming grey wolves with coyotes (<2%) and <1% to all other canid groups.The unsupervised cluster analysis was further supported by western Great Lakes population grey wolves having the lowest genetic differentiation estimates with eastern wolves (F ST = 0.06 and weighted F ST = 0.08) and western coyotes (F ST = 0.09 and 0.12), in contrast to the estimates between northern Rocky Mountains population grey wolves and eastern wolves (F ST = 0.10 and 0.10) or western coyotes (F ST = 0.12 and 0.15).

| DISCUSS ION
An estimate of the effective population size provides a means by which conservation practitioners can accurately use theory to predict forward-in-time outcomes for various viability scenarios for an endangered species (Lacy, 1995).These estimates permit one to estimate the number of generations until gene flow is required to boost the genetic diversity and concomitantly reduce inbreeding coefficients.The application of this theory to wild endangered or threatened populations has remained challenging but is centrally needed for conservation planning and simulating evolutionary outcomes (Frankham et al., 2019).One complication in the interpretation of effective population sizes is the sensitivity of these estimates to population structure (Ellegren & Galtier, 2016).Grey wolves inhabiting North America represent a diversity of demographic histories and contemporary dynamics that manifest as distinct genomic signatures.Local adaptation, compounded with social structure of grey wolves, generates population structure and increases the rate at which random genetic drift depletes their genomic variation and evolutionary potential.When geographic regions experience local extinctions from over-exploitation, dispersals will re-populate the new vacancy and genetically homogenize across proximal subpopulations over time (Ausband & Waits, 2020).Despite these recent demographic events of reintroduction or re-population, observed heterozygosity is lower than expected with significant genetic structure across the continent.
As per theory, this suggests that the effective population sizes calculated here for each grey wolf population are impacted (Ellegren & Galtier, 2016).
The comparison of the census and effective population sizes provides a more valuable metric beyond census size alone.For species with social organization, substructure, and non-random breeding, theory expects that effective population size will be a fraction of the census size (Ellegren & Galtier, 2016;Frankham, 1995).Although there are many field methods for estimating the ratio of census size to N e , these are often challenging and require an immense effort in the field.For example, using wolf Overall, the census and effective population sizes differ by approximately an order of magnitude.
We conducted a population genome-level survey of three genetic groups of grey wolves across North America and resolved deeper fine-scale resolution that was reflective of geography and demographic history.These groups correspond to the Great Lakes region, northern Quebec, and the western region of Canada and the United States.While all the populations we studied have a history of over-exploitation, each group has unique aspects to their population histories.The grey wolves of the Great Lakes carry a genetic signature of historic admixture (Heppenheimer et al., 2018;Koblmüller et al., 2009;Leonard & Wayne, 2008;Rutledge et al., 2010;vonHoldt et al., 2011vonHoldt et al., , 2016)), and habitat loss has been of consequence to wolves in northern Quebec (Larivière et al., 2000).The genetic cluster composed of the continent's western region is likely due to the shared ancestry when wolves were translocated from west-central Canada as founders for the populations in the northern Rocky Mountains with recent dispersal across the region (Hendricks, Schweizer, Harrigan, et al., 2019;vonHoldt et al., 2010).

| Northern Rocky Mountain grey wolves have declining genetic diversity
Grey wolves were restored in the northern Rocky Mountains the first evaluation of genetic structure, diversity, and connectivity over the initial 10-year recovery period (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004) inferred from microsatellite markers and reported no immediate concerns for genetic variability.However, genome sequencing advances have provided the grey wolf with a plethora of new genetic methods that avoid some central and limiting concerns when using microsatellite markers (Väli et al., 2008).As such, we encourage genetic surveys of grey wolves to consider a genome-wide reduced representation or targeted sequence-based method for large-scale population studies, which is feasible for any sample type and is less prone to calibration and ascertainment concerns of microsatellites collected across facilities, platforms, and research groups (Bonin et al., 2004;Pompanon et al., 2005).
We found genetic evidence of dispersal patterns in the Pacific Northwest, where genetic signatures clearly identified that these western continental wolf populations relied upon male-mediated dispersal for gene flow.We also detected signatures that female wolves across the western USA and southwestern Canada were significantly more differentiated from each other than males.In contrast, this pattern was not found in the females of the Great Lakes and Ontario region, likely an interaction between the population never being fully eradicated and an evolutionary history of genetic admixture with coyotes.Further, we report evidence of both significantly lower levels of genomic diversity in the northern Rocky Mountains paired with eroding diversity and higher inbreeding coefficients since 1990, explained in part by our new effective population size estimates.This temporal decline in genetic diversity was not found in the western Great Lakes wolves.
One limitation is that our genetic focus does not explore the fitness effects of such trends; however, such metrics are often central in conservation strategies.Although we currently do not report on fitness-related consequences, evaluations of such have been conducted on highly bottlenecked and inbred populations like Isle Royale and Scandinavia (Åkesson et al., 2022;Hagenblad et al., 2009;Robinson et al., 2019).The wolves of the northern Rocky Mountains currently have an increased mortality rate due to relaxed regulation.Notwithstanding, grey wolf life history of short time to sexual maturity, large litters, and dispersal can mitigate population-level risks from human-related mortality (Adams et al., 2008;Fuller et al., 2003).However, Cassidy et al. (2022) recently found significant effects of human-caused mortality on other important biological processes in wolves (e.g.pack persistence and pup production) that have implications for breeding

| Great Lakes grey wolves have a unique demographic history
Following theoretical expectations, the level of genetic richness and uniqueness is correlated with the western Great Lakes wolf demographic history of colonization and admixture (Allendorf et al., 2001).In agreement with previous findings, western Great Lakes wolves carry the lowest levels of inbreeding and the highest levels of allelic richness and private alleles.This is explained by their historic genetic exchange with other sympatric canid lineages, supported by both genetic cluster analysis and the lowest genetic differentiation with eastern wolves (F ST = 0.06 and weighted F ST = 0.08) and western coyotes (F ST = 0.09 and 0.12), in contrast to the estimates between northern Rocky Mountains population grey wolves and eastern wolves (F ST = 0.10 and 0.10) or western coyotes (F ST = 0.12 and 0.15).This demography is unique and provides an immediate mechanism by which these populations can respond to a rapidly changing world in terms of both climate and anthropogenic activity (Carmichael et al., 2008;Kagawa & Seehausen, 2020;Ottenburghs, 2021;Pacheco et al., 2022;Rius & Darling, 2014;vonHoldt et al., 2017).

| Conservation decisions in light of effective size estimates
We compiled reported population sizes across the states that compose the northern Rocky Mountains and western Great Lakes population between 1982 and 2015 from public data and found that grey wolf effective population sizes were 5.2-9.3% of the census size.Peterson et al. (1998) used demographic models of N e for Isle Royale and estimated an N e /N ratio of 16%.Further, many wild canid species will avoid mating with relatives (Ausband, 2022;Geffen et al., 2011;Sparkman et al., 2012;vonHoldt et al., 2008), and this inbreeding avoidance mechanism will increase N e .Our estimates are comparable to those for the cooperative breeding African wild dog (Lycaon pictus) where effective population sizes are 8.7-11.3% of the census size (Marsden et al., 2012).According to international conservation goals of the '50/500 rule', the genetic consequences of population subdivision are strongest in small (N e < 500) isolated populations where inbreeding depression occurs, and genomic diversity erodes due to drift.Thus, successful short-term conservation efforts can target N e ~ 50 but should target N e > 500 for the long-term survival of a species (Caballero et al., 2017;Frankham et al., 2014;Jamieson & Allendorf, 2012;Pérez-Pereira et al., 2022).As per this rule, we show that grey wolves fall above minimum effective population sizes needed to avoid extinction due to inbreeding depression in the short-term but face long-term risk of extinction on their own given their present-day effective population sizes (N e ~ 142.7-226.3).
A similar situation was also found for Scandinavian wolves, with real- ized N e below advised conservation goals (Laikre et al., 2016).Their  et al., 2020;Hendricks, Schweizer, & Wayne, 2019;Schweizer et al., 2016).The suggested effective migrant strategy would require more consideration of regional signatures of adaptive variation (Carroll et al., 2020).We envision this study as a baseline for future assessments.

| Genetic conservation of grey wolves
Species recovery plans are constructed around a core conservation biology framework referred to as 'The Three R's' (representation, resiliency, and redundancy) for reducing the risk of extinction (Shaffer & Stein, 2000).Under the ESA, this can be satisfied by maintaining multiple large, genetically robust populations across the historic range that are self-sustaining.Grey wolves have already met many of these aspects, with several populations found across the United States, and natural dispersal occurring to help occupy portions of their historic range, although the species still only occupies approximately 10-15% of its historical range (Carroll et al., 2006).With fluctuating federal protection, populations can recover, be delisted, experience reductions through human-caused mortality, and then return to federal protection, thus restarting the cycle.In addition to jurisdictional issues within the United States (Smith et al., 2016), there are also international challenges.Both populations considered here are part of a larger grey wolf population that is distributed across the United States and Canada border, making their conservation status dependent upon biological and social conditions in both countries.Joint USA-Canada conservation plans and actions have been successfully executed in the past (Bangs & Fritts, 1996), but international coordination can be complicated to maintain (Quevedo et al., 2019).Any disruption of dispersal across this international line, or decline in one country, would impact the population viability of All authors contributed towards the preparation and writing of the manuscript.

ACK N O WLE D G E M ENTS
We acknowledge the following for samples used in this study: the Mexican wolf C. l. baileyi subspecies) regions.Grey wolves in the northern Rocky Mountains were extirpated by the 1920s and were listed under the Endangered Species Act (ESA) in 1973.As such, all grey wolves in the lower 48 United States range were listed as endangered, with the exception of grey wolves living in Minnesota that were listed as threatened.The northern Rocky Mountain Wolf Recovery Plan (NRMWRP) outlined grey wolf recovery by supporting natural colonization and translocation of 66 wolves from Alberta and British Columbia to central Idaho and Wyoming's YNP during the winters of 1995 and 1996 (59 FR 60266; U.S. Fish and Wildlife Service, 1987).Dispersers from YNP expanded into adjacent Montana, Idaho, and Wyoming counties (collectively referred to as the Greater Yellowstone Ecosystem), and dispersers from central Idaho expanded into adjacent Montana, Wyoming, and Oregon.Beginning in the late 1990s, periodic dispersing wolves from southern British Columbia and the northern Rocky Mountains were documented in the Pacific Northwest states of Washington, Oregon, and northern California.By 2011, the first wolf entered Oregon with confirmed reproduction in 2015.The western Great Lakes population is composed of the eastern portion of the Dakotas, Minnesota, Iowa, Wisconsin, a northern portion of Illinois, and Michigan (lower and upper peninsula).Grey wolves in Minnesota were first protected under the ESA in 1974, with subsequent expansion into Wisconsin and Michigan by the early 1990s genome ) of2,202,059,258 and 123,842,264  nucleotides for autosomes and the X chromosome, respectively.The geographic region was used as an identifier for the function homozyg in PLINK v1.9 ously published canids: 106 reference western coyotes (C.latrans) from vonHoldt et al. (2022) and 30 reference eastern wolves (C.lycaon) from Heppenheimer et al. (2018) (Table and repeated the filtering.The result is a dataset of 162,109 minimally filtered SNPs genotyped in 437 grey wolves from Canada (n = 92) and the United States (n = 345), with a subset of 81,595 loci referred to as the 'statistically neutral and unlinked' SNPs.We inferred sex for individuals bioinformatically based on the depth of reads mapped to the Y chromosome.Of the 437 wolves, field-based observations identified 104 females and 118 males.

Figure
Figure S2).Three genetic clusters broadly represent three distinct demographic histories of western, central, and eastern regions of the continent.We divided the western cluster into two subclusters, one to reflect the shared demography of southwestern Canada and western USA through the translocation and colonization of wolves in the northern Rocky Mountains population, and the other representing northern Canada (TableS1).The other two clusters represent north-

Finer
-scale clustering revealed a stronger role of geographic isolation, with more resolution of substructure within USA's northern Rocky Mountains and the Pacific Northwest regions (Figure 1c).The shared assignments across three genetic partitions reflect the shared genetic ancestry across large geographic distances due to the translocation of grey wolves in 1995 and 1996 (British Columbia, Alberta, and Montana) to central Idaho and the Greater Yellowstone Ecosystem (mean Q: partition 1 BC = 0.43, ID = 0.14, GYE = 0.22; partition 2 BC = 0.25, ID = 0.40, GYE = 0.07; partition 3 BC = 0.09, ID = 0.13, GYE = 0.65).Populations with shared demographic histories (northern Canada vs. western USA and southwestern Canada, F ST = 0.034) levels found among northern/southwestern Canada and the western USA regions (H O = 0.223 and 0.220), and the Great Lakes TA B L E 1 Average and weighted Weir and Cockerham estimates (above and below diagonal, respectively) of genetic differentiation (F ST ) across 81,595 SNPs between geographic regions of grey wolves (see Figure 1a for population abbreviations).
the northern Rocky Mountains (n = 188) and the western Great Lakes(n = 199).The preceding analysis identified the distinctiveness between the northern Rocky Mountains and western Great Lakes population segments as per their divergent assignment probabilities (K = 3 and K = 9) (Figure1b,c, FigureS2).We found that six (4.5%) of the northern Rocky Mountains wolves had assignments to a cluster divergent from their geographic origins at K = 3 (when Q > 0.00001, Q = 0.01-0.25),all of which were individuals sampled in the Pacific Northwest.The misclassification of western Great Lakes wolves is more varied due to assignments to the proximate Canada wolf populations at K = 3 (Q = 0.01-0.86).This pattern continued at K = 9, where the highest non-Rocky Mountains assignments were wolves assigned to Canada's Northwest Territories Province (Q = 0.01-0.37),concordant with a shared demographic history.We identified seven western Great Lakes individuals with assignments (several samples in Isle Royale NP, Q = 0.01-0

western
Great Lakes, Michigan (including wolves on Isle Royale) had the lowest estimates (H O = 0.219) compared to Minnesota (0.225) and Wisconsin (0.231).We restricted the analysis to samples only with known years of sample collection between 1990 and 2020 within the population of the northern Rocky Mountains (n = 137) and western Great Lakes (n = 86) to survey changes in diversity over time.Using Pearson's product-moment correlation, we found that all heterozygosity estimates for the northern Rocky Mountains population significantly declined over the 30 years surveyed (H O : R = −.41,p = 8.3 × 10 −7 ; tions (approximately 200 years) from autosomal SNPs for each of F I G U R E 2 Heterozygosity (observed and expected) trends for the (a) northern Rocky Mountain (n = 137) and (b) western Great Lakes (n = 86) distinct population segments in the United States for a 30-year period between 1990 and 2020 (Y-axis).Pearson correlation coefficients and significance values are provided.Shaded area indicates the 95% confidence interval.
5 wolves 36 generations (144 years) ago, with western USA/southwestern Canada next largest for estimates of 1989.4 wolves 41 generations (164 years) ago, then Great Lakes/Ontario with 878.7 wolves (45 generations or 180 years ago), and finally northern Quebec at low estimates maxing at 464.8 wolves 46 generations (184 years) ago.
displayed a steep and rapid effective rate of loss (m = −45.6)per generation while the western Great Lakes population's decline was shallower (m = −14.4)(Figure 3a).The northern Rocky Mountains experienced a dramatic shift 20 generations ago losing 72.8 wolves per generation.In that same time frame, the western Great Lakes was losing 4.0 wolves per generation.Their current-day respective estimates are N e_RM = 141.7 and N e_GL = 226.3,after having effective population size estimates reduced by 1928.6 and 542.1 Both regional populations have a history of substantial expansion in census population sizes between 1982 and 2010 when the northern Rocky Mountains were estimated to have N ~ 1723 and western Great Lakes at N ~ 4321 wolves, remaining mostly stable to the present-day estimates of N ~ 1881 and 3025, respectively (Figure 3b).We estimated that the western Great Lakes effective population size has remained stable since 1990 with an average rate of growth larger than that of the northern Rocky Mountains (GL m = 0.21; RM m = −0.05),with significantly higher effective population size estimates for western Great Lakes (N e = 226.6)than the northern Rocky Mountains (N e = 143.8)(t-test unequal variance p = 1.420 × 10 −11 ).Lastly, we estimated the temporal trend of N e /N collectively for the northern Rocky Mountains and the western Great Lakes and found the effective population size remained at 5.2-9.3% of the census size since mid-2000s (Figure 3b).We estimated that the decay in heterozygosity for the northern Rocky Mountains had an initial level of H O ~ 0.235 in 1991 and decayed to 0.208 by 2020 (approximately eight generations) genotyped in 465 canids (179 grey wolves from the northern Rocky F I G U R E 3 Locally estimated scatterplot smoothed (loess) trend lines of population effective size (N e ) histories for (a) each of the four identified regional genetic clusters and the regional populations in the United States.The vertical dashed line in each panel indicates the acceptance of the U.S. Endangered Species Act into law in 1973.(b) Observed (N) and inferred population effective size (N e ) histories for the northern Rocky Mountain and the western Great Lakes populations in the United States.We assumed 4 years per generation.The inset displays the ratio of N e to N since 1982-2015 for each of the two populations with values included.
dispersal and density data on the Perch Lake pack (N m = 5, N f = 5) in Minnesota, Chepko-Sade et al. (1987) estimated effective population size with two methods: the root mean square (variance) method (N e = 804) and the 85th percentile distance of the original dispersal distribution method (N e = 1660.7).In comparison, we provided a genomic N e estimate of 222.6 wolves in 1987 for the western Great Lakes, roughly 13-28% of that derived from wolf dispersal and density data.Further, earlier population estimates from 26 microsatellite data of Yellowstone National Park wolves reported N e ranging between 6 and 22.6 for 1995-2004 and the respective census sizes of 21 and 80 (range N e /N = 0.10-0.37)(vonHoldt et al., 2008).Genomic-based inferences still face challenges albeit different from field-based inferences; regardless, estimates are critical for shaping appropriate conservation management plans.Understanding this relationship is important because management applies to actual populations which are observed and managed based on census size, not effective population size.Using genomic data from these populations, we show that this ratio is different in different parts of the distribution.
through a reintroduction programme in the mid-1990s and a handful of dispersing wolves southward from Canada into northwestern Montana, which successfully established several that contributed towards the first of many delisting proposals for this population in 2003.A study by vonHoldt et al. (2010) provided and gene flow.Given the difficulty states have faced in meeting their goals of significant population reduction (e.g.Idaho's goal of 500 wolves with an estimated 1270 census size, Idaho Fish and Game Grey wolf management plan draft January 2023), the effective population size estimates are then interpreted to be strongly influenced by the number of breeding wolves and gene flow, less from census size.Current management actions that seek to reduce overall populations and permit hunting during the breeding season have the greatest potential to have negative consequences on effective population sizes.
ultimate suggestion was to increase N e and promote methods that would increase genetic exchange via 3-5 effective migrants per generation with neighbouring populations.Notably, such goals are clearly possible within the ESA framework which defines 'conservation' in section 3 to include 'the use of all methods and procedures which are necessary to bring any endangered species or threatened species to the point at which the measures provided pursuant to this Act are no longer necessary'.There are known dispersers, albeit unknown if they are effective dispersers, between southwestern Canada and the U.S. Rocky Mountains.Combined with the shared ancestry due to translocation from the western Canada and northern Rocky Mountain grey wolf populations, demography is a core feature that shapes conservation-relevant metrics.Further, wolves in North America can originate from dramatically different regions with distinct collections of local adaptations and ecotypes (Carroll the wolves.The Assistant Secretary of the Interior is quoted, regarding the ESA that '…it is in the best interests of mankind to minimize the losses of genetic variations.The reason is simple: they are potential resources.They are keys to puzzles which we cannot solve, and may provide answers to questions we have not yet learned to ask' (H.R.Rep.No. 93-412, pp.4-5, 1973).Such Congressional intent clearly displays the intent of including all means for the conservation of genetic variation.Further, human activity homogenizes the landscape on which endangered species rely, and such activities '… threaten their -and our own -genetic heritage.The value of this genetic heritage is, quite literately, incalculable' (93D Congress Report, 1st Session, No 93-412, page 143).The minimum effective population size of 500 necessary to ensure long-term population viability has been difficult to apply in practice.There are many reasons for this.One reason is the abstractness-it can be hard for a manager to know what the effective population size of the population they are managing is when what they can count is the census size.In 2021, the northern Rocky Mountains had a census size estimated at 3354 and western Great Lakes at 4526.However, we can then translate these values to an effective population size ranging between 201 and 335 wolves for the northern Rocky Mountains and 272 and 453 for the western Great Lakes.Given the strong skew in the effective-to-census size ratio in grey wolves, larger wolf populations are necessary to ensure long-term adaptation and survival.Disperser success is an additional critical factor for long-term survival of the species, promoting gene flow that will reduce inbreeding and elevate effective population sizes through increased allelic variation and demographic rescue (Newmark et al., 2023).Dispersers are often challenged by utilizing lower quality corridors with high mortality risk to find suitable areas for establishing new territories (Oakleaf et al., 2010).The protection of grey wolf dispersers between wolf populations is thus important to improve their effective population sizes for long-term persistence.AUTH O R CO NTR I B UTI O N S Bridgett M. vonHoldt, Daniel R. Stahler, and Robert K. Wayne conceived and designed the research study.Daniel R. Stahler, Marco Musiani, Rolf Peterson, John Stephenson, Kent Laudon, and Erin Meredith collected the samples.Bridgett M. vonHoldt generated the RADseq SNP data and analysed the data with guidance from Daniel R. Stahler, Kristin E. Brzeski, Rolf Peterson, and Kent Laudon.

L
. David Mech, Michel Crête, and the Michigan Department of Natural Resources.We are appreciative of comments provided by Stephen Gaughran, Armando Caballero, and all reviewers who greatly improved our manuscript.We are grateful to the Turner Endangered Species Fund for partial funding support.Daniel Stahler and Yellowstone Wolf Project efforts were funded by NSF DEB-0613730 and DEB-1245373, Yellowstone Forever, and their many donors, especially Valeria Gates and Annie and Bob Graham.This research was also implemented under the National Recovery and Resilience Plan (NRRP), Project title 'National Biodiversity Future Center -NBFC'.CUP J33C22001190001.
.56) to Canada's Northwest Territories Province, two assigned to Idaho (sampled in MN and WI, Q = 0.99), Although we found that the northern Rocky Mountains and western Great Lakes populations carried comparable observed heterozy-gosity levels (H O , H E =[0.211, 0.224] and [0.211, 0.211], respectively), the per-state composition was quite variable (Table2b).Estimations at the state level revealed that in the northern Rocky Mountains, the four samples from California were the most genetically diverse (H O = 0.562), followed by Montana (0.333), Washington (0.298), Oregon (0.285), Idaho (0.245), and Wyoming (0.238) (Table2b).In the TA B L E 2 Average expected and observed heterozygosity (H E and H O , respectively) and effective population size (N e from past 50 generations) estimates for each (a) major geographic location (p-values are from a Welch two-sample t-test of unequal variance between H E and H O ) and (b) regional population within the United States.Diversity estimates were derived from the statistically neutral SNP set while effective population size estimates from the minimally filtered SNP set.Abbreviation: n, sample size.a Includes grey wolves from Isle Royale National Park in Lake Superior.