A Dominant Founder Lineage Has Possible Fitness Costs for the Endangered Mexican Grey Wolf

TL;DR

The Mexican grey wolf population faces genetic challenges due to inbreeding, but interbreeding strategies have helped reduce these issues and improve genetic diversity.

Contribution

The study provides genomic insights into inbreeding and lineage contributions in the Mexican grey wolf ex situ breeding program.

Findings

Inbreeding levels decreased after inter-lineage pairings, possibly due to fitness costs of inbreeding depression.

Higher inbreeding and autozygosity were linked to reduced lifespan and reproductive performance.

Ghost Ranch and Aragón lineages helped reduce inbreeding in the genetically distinct McBride lineage.

Abstract

The Mexican grey wolf (
Canis lupus baileyi
) is an endangered and genetically distinct subspecies of grey wolf adapted to the warm climates of the southwestern United States and Mexico. Following centuries of eradication efforts, Mexican grey wolves were protected under the Endangered Species Act in 1976, prompting an international ex situ breeding program to preserve their genetic legacy. Seven remnant wolves founded the three distinct lineages of this program: McBride, Ghost Ranch, and Aragón. However, concerns of rising inbreeding levels motivated the implementation of inter‐lineage pairings to decelerate genetic erosion. To evaluate the genetic health of the ex situ managed population, we analysed genome‐wide variation across several generations of breeding wolves and their descendants (n = 179). Despite ongoing conservation interventions, we found a sharp decline in effective population size over the past 50 generations. However, estimates of inbreeding diminished over time following the inter‐lineage pairings, reflecting efficient pedigree management or potential fitness costs of inbreeding depression, as higher levels of inbreeding and autozygosity were associated with reduced lifespan and reproductive performance. Increasing contributions of Ghost Ranch and Aragón ancestry also aided in reducing levels of inbreeding observed in the more genetically distinct McBride lineage. These findings highlight the challenges of maintaining genetic diversity in bottlenecked populations, reveal fine‐scale patterns of inbreeding within the pedigree and underscore the urgent need to preserve genetic contributions from all founder lineages. Genomic insights from this study can guide breeding strategies to reduce kinship, balance ancestry and mitigate inbreeding depression, safeguarding the evolutionary legacy and the long‐term persistence of the Mexican grey wolf.

Figures

(A) Pedigree of the Mexican grey wolf studbook, centred on individuals sequenced in this study (visualized using R's kinship package, Sinnwell et al. 2014; shaded symbols; squares = males; circles = females; double horizontal lines = breeding between first‐degree relatives; dashed lines = links between sub‐pedigrees; studbook ID shown below each symbol). (B) Estimated effective population size (N e) over the past 50 generations (4 years/generation) with LOESS‐smoothed trend and 95% confidence interval. Vertical dashed lines mark two major events (ESA = Endangered Species Act). (C) Yearly mean and standard error estimates of inbreeding and autozygosity within the ex situ breeding population with a linear model trendline for pre and post‐ 1995 inter‐lineage crossed populations (vertical dashed line) provided. Relationships as a linear regression with trend lines (with shaded 95% confidence intervals) between the number of homozygosity segments (ROH) and (D) lifespan (age at death) and (E) the number of pups, and (F) the average ROH length as a function of the number of pups produced by an individual.

(A) PCA of unlinked SNPs, with points coloured by maximum inferred ancestry assignment to the founder lineages. (B) Inbreeding coefficient (F ROH), (C) distribution of run‐of‐homozygosity lengths (Kb) and (D) decay of linkage disequilibrium for each founder lineage assignment (standard error depicted as the shading around the mean line).

(A) Triangle plot of the supervised admixture assignment probability (Q) for each query individual with respect to the three reference founder lineages (posterior probability lineage assignment indicated by colour). (B) The number of breeding individuals assigned to each lineage ancestry and the (C) mean measures of autozygosity (F ROH) for each assigned lineage ancestry within the dataset across each year (vertical dashed line marks the 1995 inter‐lineage crossing start date). (D) Average number of heterozygous genotypes with predicted deleterious variants across each assigned lineage ancestry.

Tables

TABLE 1. Timeline for the Mexican grey wolf.

Time	Event	References
Early/mid‐1900s	Widespread predator eradication programs began in the US	Mech and Boitani (2003)
1959–1961	‘Ghost Ranch’ Mexican wolf ex situ lineage established in New Mexico	Hedrick et al. (1997)
1965	‘Aragón’ Mexican wolf ex situ lineage established in Mexico	Hedrick et al. (1997)
1970	Mexican wolf extirpation from the United States	Brown (1983)
1976	Mexican wolf listed as endangered under the Endangered Species Act	Hedrick et al. (1997)
1979	The original Mexican Wolf Recovery Team was formed with representatives from Mexico, Arizona, New Mexico and the FWS	Harding et al. (2016)
1977–1980	‘McBride’ ex situ lineage established in captivity from animals captured in the wild in Mexico to preserve remaining genetic diversity	Hedrick et al. (1997); García‐Moreno et al. (1996)
1980	Mexican wolf extirpation from Mexico	García‐Moreno et al. (1996)
1982	FWS wrote the first iteration of the Mexican Wolf Recovery Plan, which was then signed by representatives from the United States and Mexico.	U.S. Fish and Wildlife Service (1982)
1990	FWS began collaboration with Museum of Southwestern Biology to archive Mexican wolf Recovery Plan mortalities and serial blood samples	Museum of Southwestern Biology (n.d.)
1995	Inter‐lineage crossing of founder lineages	Hedrick and Fredrickson (2008)
1998	The first 11 ex situ‐raised Mexican wolves were released into the Blue Range wolf recovery area of Arizona	U.S. Fish and Wildlife Service (1998), Harding et al. (2016)
2011	Subsequent collaboration with Mexico resulted in the release of ex situ wolves in northern Mexico	Harding et al. (2016)
2015	Mexican wolves were formally recognised as a distinct subspecies of grey wolf ( Canis lupus baileyi )	Fredrickson et al. (2015)
2015	Rule to expand recovery area for the experimental population.	US Fish and Wildlife Service (2015)
2017	Finalisation of Mexican Wolf Recovery Plan for long‐term recovery	US Fish and Wildlife Service (2017)
2022	The FWS published a second revision of the Mexican Wolf Recovery Plan	US Fish and Wildlife Service (2022)
2024, 2025	Updates on the progress and challenges posed by the Mexican Wolf Recovery Plan emphasising the need for genomic estimates of diversity	US Fish and Wildlife Service (2024, 2025)

TABLE 2. Counts (proportions out of total ROH blocks per group in parenthesis) for each ROH block size category per founder ancestry or lineage group.

Group	1–2 Mb	2–4 Mb	4–8 Mb	8–16 Mb	> 16 Mb
Aragón (n = 54)	151 (0.08)	410 (0.20)	630 (0.31)	530 (0.26)	297 (0.15)
Ghost Ranch (n = 22)	60 (0.09)	129 (0.20)	207 (0.32)	143 (0.22)	117 (0.18)
McBride (n = 103)	423 (0.08)	826 (0.16)	1521 (0.29)	1444 (0.27)	1065 (0.20)
Founders (n = 27)	97 (0.07)	183 (0.14)	394 (0.29)	353 (0.26)	325 (0.24)
Inter‐lineage crosses (n = 152)	537 (0.08)	1182 (0.18)	1964 (0.30)	1764 (0.27)	1154 (0.18)

Full text

Introduction

1

Preserving locally adapted evolutionary lineages is critical for maintaining functional biodiversity. Such lineages often warrant formal subspecies taxonomic designation (Dufresnes et al. 2023) to recognise the valuable subset of unique phenotypic and genomic variation while still carrying the characteristics that distinguish the species (Ryder 1986; Phillimore and Owens 2006; Coates et al. 2018). These subspecies represent one of several building blocks included in species conservation plans, guided by the larger conservation goal of preserving adaptive variation and diversity within polytypic species. This design is formally recognised as ‘The Three R's’ of conservation, serving the fundamental mission of the Endangered Species Act (ESA) of 1973 to maintain: (1) resilience for the species to withstand and recover from disruption or threats; (2) redundancy for the species to exist in multiple, independent entities; and (3) representation, or the ability of a species to adapt by valuing the full spectrum of variation (Shaffer and Stein 2000). These values emphasise the crucial role that both abundance and adaptive capacity play in the long‐term recovery of species of conservation concern.

When distinct evolutionary lineages face decimation, their decline in abundance may directly impact life history and demographic aspects of their viability (Wright 1931; Hohenlohe et al. 2020; vonHoldt et al. 2025). This interaction is amplified for social species with reduced effective population sizes (N e) from skewed reproduction and large ranges, exacerbating threats to genetic diversity and adaptive potential (Wright 1931; Frankham 1995; Lynch and Lande 1998; Frankham et al. 2002; Charlesworth 2009; Hohenlohe et al. 2020; vonHoldt et al. 2025). Following population genetic theory, fragmented populations that experience severe isolation, loss of abundance and inbreeding can rapidly escalate the risk of extinction due to an increased genetic load of alleles with deleterious impact, accelerating the loss of individual‐level fitness known as inbreeding depression (Lande 1994; Hedrick and Kalinowski 2000; Keller and Waller 2002; Chapman et al. 2009; Charlesworth and Willis 2009; Frankham et al. 2014).

Faced with the risk of extinction and insufficient time to address the underlying causes in the wild, conservationists may translocate individuals of a species, subspecies or distinct evolutionary lineage to ex situ facilities to help preserve viability and adaptive genetic variation via captive breeding programs (Canessa et al. 2015). Nevertheless, the challenge persists as this newly bottlenecked population established with a handful of founders will continue to face the same major threats to small, isolated endangered populations (Lynch and Lande 1998). Conservationists attempt to address these obstacles by managing ex situ populations using a studbook that tracks pedigree relationships and guides breeding to reduce kinship and preserve genetic diversity (Lacy 1995, 2012; Lacy et al. 1995; Fernández et al. 2001, 2003, 2004; Pelletier et al. 2009; Ballou et al. 2010; Willis and Willis 2019; Galla et al. 2020, 2021).

While many wild animals, particularly social species, naturally display inbreeding avoidance, such mate choice options are not always permitted or possible when under human care for the survival of the species (Pusey and Wolf 1996; Penn and Potts 1999; Tregenza and Wedell 2001; Geffen et al. 2011; Szulkin et al. 2013; Duthie and Reid 2016). Expectations from extensive research suggest the most significant consequences of this inbreeding are found in reduced survival and fertility of offspring produced from inbreeding events (Charlesworth and Charlesworth 1987; Pusey and Wolf 1996; Hedrick and Kalinowski 2000; Keller and Waller 2002; Wright et al. 2008; Charlesworth and Willis 2009; Huisman 2017; Kardos et al. 2016; Sumreddee et al. 2020). However, studies have also shown that these fitness effects can be genetically rescued by gene flow or strategic intercrossing of inbred lineages (Coulson et al. 1998; Hogg et al. 2006; Lippman and Zamir 2007; Frankham 2015; Whiteley et al. 2015; Hedrick and García‐Dorado 2016; Hasselgren et al. 2018). Thus, it is especially critical for the ex situ management of endangered species and distinct evolutionary lineages to understand the degree to which the loss of genetic variation, small N e and population isolation synergistically impact individual‐level fitness and functional diversity.

The use of genetic tools can offer deeper insights into the genetic stability of an ex situ population, as undocumented breeding events, unsampled pedigree contributors or inbreeding can modify the level of shared genetic variation expected solely from a studbook pedigree (Russello and Amato 2004; Witzenberger and Hochkirch 2011; Jiménez‐Mena et al. 2016; Galla et al. 2020, 2021; Wright et al. 2021). For example, the pedigree‐based inbreeding coefficient (F) is the probability that two alleles at a locus are identical by descent (IBD), often defined relative to the pedigree founders that are assumed to be unrelated (Wright 1931; Powell et al. 2010). However, pedigree‐based F estimates may be underestimated if calculated from studbook information that fails to capture signatures of undocumented inbreeding events or other unobserved factors, such as random Mendelian segregation and founder relatedness (Oliehoek and Bijma 2009; Pelletier et al. 2009; Knief et al. 2015; Wang 2016). Therefore, there is a growing consensus to incorporate complementary genome‐wide estimates of inbreeding derived from comparisons between observed levels of homozygosity in the population to what is neutrally expected (F), as well as from characterising the autozygous segments of the genome (F ROH) (Wright 1931; Powell et al. 2010; McQuillan et al. 2008). These realised estimates can better inform management decisions regarding inbreeding avoidance and the retention of genetic diversity.

Here, we explored the genomic profiles and viability of the ex situ managed population of Mexican grey wolf ( Canis lupus baileyi ), an endangered subspecies endemic to North America. Uniquely adapted to the southwestern United States and northern and central Mexico, this subspecies was nearly driven to extinction by centuries of habitat loss, hunting and federally sponsored predator eradication efforts (Parsons 1996, McBride 1980, Brown 1983) (Tables 1, S1). After receiving protection under the Endangered Species Act in 1976 (USFWS 1976), seven remnant Mexican grey wolves from zoos (Ghost Ranch n = 2 and Aragón n = 2) or wild translocations (McBride n = 3) established the three founder lineages of a bi‐national captive ex situ breeding program designed to preserve the genetic legacy of this distinct evolutionary lineage (Hedrick et al. 1997). Advances in genomic technology confirmed the Mexican grey wolf as a distinct subspecies of grey wolf with minimal introgression from other canids (National Academies of Sciences, Engineering, and Medicine 2019 and references within), yet there were no other confirmed sightings of Mexican grey wolves remaining in the wild.

Extirpated from the wild with no hope for the discovery of new founder lineages, preserving genetic diversity within the ex situ population became a critical focus of the recovery program. To mitigate the effects of inbreeding depression within the small founder populations, inter‐lineage crossbreeding was initiated in 1995 and maintained through a cross‐institutional studbook (Fredrickson and Hedrick 2002). By 1998, wolves with high mean kinship and strong representation in the ex situ population were used to establish an experimental population in the wild (Hedrick and Fredrickson 2008; Hendricks et al. 2016). Nevertheless, comprehensive genomic studies remain essential to quantitatively evaluate the preservation of genetic diversity and extent of inbreeding in the ex situ breeding and reintroduced populations. Early genetic studies have utilised variation in allozyme loci and mtDNA sequences to confirm the genetic distinctiveness and subspecies identity of Mexican grey wolves (Shields et al. 1987; Wayne et al. 1992), as well as multilocus DNA fingerprinting to assess the evolutionary placement and genetic diversity of the Mexican grey wolf subspecies and founder lineages (Hedrick 1971; Wayne 1995; García‐Moreno et al. 1996; Hedrick et al. 1997). Yet, despite the decreasing costs of genetic surveillance, there have been no efforts to characterise genome‐wide patterns of inbreeding for this ex situ population, relying on pedigree‐based estimates of inbreeding for the management of Mexican wolf founder lineages (Kalinowski et al. 1999; Fredrickson and Hedrick 2002; Hedrick and Fredrickson 2008; Clement et al. 2024).

To evaluate the genetic viability of the ex situ population following generations of managed captive breeding among genetically bottlenecked founders, we analysed reduced‐representation genome‐wide genotype data from a subset of the pedigreed breeding individuals archived in a museum biorepository that temporally spans the pedigree before and after the inter‐lineage crosses. Although reproductive outcomes rely on pedigree‐based management decisions intended to limit the loss of diversity, we aimed to provide an integrated assessment of genomic diversity, founder ancestry, inbreeding and fecundity to inform future management strategies for the genetic stewardship of this and other critically endangered organisms.

Materials and Methods

2

Samples, DNA Extraction and Restriction Site‐Associated DNA Library Construction

2.1

We received tissue samples from 187 breeding Mexican grey wolves from the Museum of Southwestern Biology, Divisions of Mammals and Genomic Resources, University of New Mexico (loan 2025.058.Mamm) that represent a subset of the entire Mexican grey wolf studbook of 2751 individuals born before 2019, along with their meta‐ and life history data for pedigree‐based comparisons, although cause of death was not included in the shared documents (Table S1; Figure 1A). We followed the manufacturer's protocol to obtain genomic DNA using the Qiagen DNEasy Blood and Tissue kit (Qiagen), which we quantified using a Qubit fluorometer system (Thermo Fisher Scientific). We prepared the genomic DNA for restriction site‐associated DNA‐sequencing (RADseq) following Ali et al. (2015). Briefly, we digested 75 ng of genomic DNA with the Sfb1‐HF restriction enzyme (New England Biolabs). We ligated a unique 8 bp biotinylated BestRAD adapter to each specimen to permit downstream pooling of up to 48 individual libraries. We used Agencourt AMPure XP Beads (Beckman Coulter) for purification and size selection steps for fragments after sonication to 300–400 bp on a LE220 Covaris at the Princeton University Lewis Sigler Genomics Core facility. We used Streptavidin Dynabeads M‐280 (Invitrogen) to enrich each genomic library for adapter‐ligated fragments that carried the biotinylated adaptors, which were then prepared using the NEBNext Ultra II DNA Library Prep Kit (New England Biolabs) to repair ends, ligate the Illumina P2 i5 adapter and amplify to incorporate a second unique barcode for each pooled library. We sequenced the prepared genomic libraries in the 2 × 150nt configuration on a NovaSeq at the Yale Center for Genome Analysis to collect 4 million reads per sample.

(A) Pedigree of the Mexican grey wolf studbook, centred on individuals sequenced in this study (visualized using R's kinship package, Sinnwell et al. 2014; shaded symbols; squares = males; circles = females; double horizontal lines = breeding between first‐degree relatives; dashed lines = links between sub‐pedigrees; studbook ID shown below each symbol). (B) Estimated effective population size (N e) over the past 50 generations (4 years/generation) with LOESS‐smoothed trend and 95% confidence interval. Vertical dashed lines mark two major events (ESA = Endangered Species Act). (C) Yearly mean and standard error estimates of inbreeding and autozygosity within the ex situ breeding population with a linear model trendline for pre and post‐ 1995 inter‐lineage crossed populations (vertical dashed line) provided. Relationships as a linear regression with trend lines (with shaded 95% confidence intervals) between the number of homozygosity segments (ROH) and (D) lifespan (age at death) and (E) the number of pups, and (F) the average ROH length as a function of the number of pups produced by an individual.

Bioinformatic Processing

2.2

We demultiplexed sequences per pool index using modules in STACKS v2.6 (Catchen et al. 2013; Rochette et al. 2019). First, we used the process_radtags module to retain reads that carried no more than 2 nucleotide mismatches to the expected barcode sequence and had a quality score ≥ 10. We removed duplicated reads with the clone_filter module before mapping these cleaned reads to the Canfam3.1 reference dog genome (GCA_2285.2) (Lindblad‐Toh et al. 2005) using bwa‐mem v0.7.19 (Li 2013), where we also included the Y chromosome (KP081776.1; Li et al. 2013). We converted the files to the BAM format with samtools v0.1.18 (Li et al. 2009). We discovered and annotated variant sites using the Marukilow model in the gstacks and populations modules in STACKS (–vt‐alpha and –gt‐alpha with p = 0.01) to identify and genotype SNP loci across all the samples sequenced. We utilised several filtering steps to retain a dataset of high‐confidence loci. First, we used VCFtools v0.1.17 (Danecek et al. 2011) to exclude loci that carried singleton or private doubleton alleles or had < 10× sequence coverage. We then excluded individuals with more than 20% missing data (–mind), removed loci with a minor allele frequency (MAF < 0.03) and required at least an 80% genotyping rate per locus (–geno) in PLINK v1.90b3i (Chang et al. 2015). This was referred to as the ‘minimally filtered dataset’.

Filtering for Population Demography

2.3

For demographic analyses that required neutral and unlinked loci, we constructed a ‘demographic dataset’ by removing loci in linkage disequilibrium (LD) with a stringent genotype correlation (r ^2^ > 0.2) in PLINK (–indep‐pairwise 50 5 0.2) and loci that significantly deviated from Hardy–Weinberg Equilibrium (–hwe, p = 0.001).

Estimating the Effective Population Size

2.4

We estimated the recent history of the effective population size (N e) over the past 200 generations using the LD‐based algorithm implemented in GONe (Santiago et al. 2020) on the demographic dataset, which is concordant with data expectations (Beichman et al. 2017). Accounting for genetic drift, we assumed a 1 cM recombination rate per 1 Mb across the genome for N e estimates in GONe. We utilised the unphased autosomal data option, did not include a MAF filter, permitted a maximum of 50,000 loci per chromosome and opted to ignore pairs of SNPs with a recombination rate over 0. We translated generations into years using 4 years per generation as the unit of time (Mech et al. 2016; vonHoldt et al. 2008) and the ‘current’ time as the year of death of the most recently sequenced individual in the pedigree. To show a biological ‘confidence interval’, we included the estimate for 2 and 6 years per generation. We averaged N e over the most recent eight generations to avoid biases from any lingering artefact in the first four generations (Novo et al. 2023). We focused exclusively on N e estimates over the last 50 generations (approximately 200 years), which is most pertinent to the demographic history of Mexican wolf decline and ex situ breeding histories. These estimates were compared to the census size of animals accounted for in the pedigree in each generation. The most recent pedigree year of death for this study was 2022.

Genetic Ancestry Assignment and Diversity

2.5

Using the demographic dataset, we used the admixture model (noadmix = 0, usepopinfo = 1) for posterior probability population assignment in STRUCTURE v2.3.4, permitting assignments up to two generations back (gensback = 2) at K = 3 (Pritchard et al. 2000). We implemented the use of both the popdata and popflag parameters with burnin = 10,000 and numreps = 20,000, along with all other flags at default settings. These probabilistic lineage assignments were used to categorise each individual into a single founder lineage. We then summarised the genetic structure within these individuals by performing a Principal Component Analysis (PCA) as a non‐model clustering method implemented in flashPCA v.1 (Abraham et al. 2017). We also calculated identity by descent probabilities (IBDp) across the selected range of missingness threshold values for autosomal SNPs in the demographic dataset using the –genome flag in PLINK. We clustered the dataset based on IBDp estimates and reported the trends using multidimensional scaling (MDS) in PLINK with the flags –cluster and –MDS‐plot 4, respectively. We statistically tested the clustering patterns for each coordinate with respect to founder lineage membership and life history traits. We then analysed genetic diversity within each founder lineage assigned by maximum ancestry by calculating observed homozygosity (H O), as well as expected homozygosity (H E) and the inbreeding coefficient (F) from these estimates under Hardy–Weinberg equilibrium in PLINK using the –het flag using the demographic dataset. We also estimated the proportion of founder lineage ancestry for each individual using a supervised maximum‐likelihood method in the program ADMIXTURE v1.3.0 using the demographic dataset (Alexander et al. 2009). We specified a value of K = 3 to indicate the three genetic reference populations (K) that corresponded to the three founder lineages (McBride, Ghost Ranch and Aragón). The reference populations were composed of individuals born in each of the three founder lineages before the inter‐lineage crossing design was initiated in 1995. Individuals born after 1995 were identified as query wolves to be assigned to one or more founder lineages. We statistically assessed differences in ancestry proportions among groups using a non‐parametric permutation‐based multivariate analysis of variance (permanova) in the R package vegan with 999 permutations based on Euclidean distance (Oksanen et al. 2025).

Genome‐Wide Estimates of Autozygosity and Linkage Disequilibrium

2.6

We calculated autozygosity as a genomic measure of inbreeding by estimating runs of homozygosity (ROH) in the minimally filtered dataset, which considers the presence and length of homozygous tracts within individuals, using PLINK's –homozyg flag. We identified tracts of homozygosity in the autosomes and grouped consensus segments across individuals, allowed ≤ 3 heterozygous calls and 50 missing calls in a window with up to 200 loci per scanning window (parameters: –dog, –chr 1–38, –homozyg group, –homozyg‐window‐het 3, –homozyg‐window‐missing 50, –homozyg‐window‐snp 200, –missing‐genotype 0; McQuillan et al. 2008; Ceballos, Hazelhurst, and Ramsay 2018). We converted ROH values to an individual inbreeding coefficient (F ROH) by dividing the summation of the total length of ROH per genome by the total length of the autosomal genome (McQuillan et al. 2008). We statistically tested patterns of F ROH across maximum founder ancestry assignments, as well as across lifespan. After ROH identification, we assessed lineage‐specific differences by grouping ROH blocks into size‐based categories: 1–2, 2–4, 4–8, 8–16 and > 16 Mb. The shorter categories represent expected sizes after a longer history of recombination and greater genomic diversity, while longer ROH sizes are expected due to more recent inbreeding events and greater overall genomic IBD (Ceballos, Joshi, et al. 2018). We also estimated linkage disequilibrium (LD) with the minimally filtered dataset using PLINK to calculate genotypic correlation as a proxy for LD using the –r^2^ flag across 1000 Kb windows within the following groups of individuals: cross‐lineage versus founder individuals; and comparing across the three founder lineages. If loci were separated by > 100 Kb, they were treated as different linkage blocks for plotting.

Associations With Fitness

2.7

Using the minimally filtered dataset, we excluded surviving individuals who had not yet reached their total reproductive output and associated various genomic metrics with two proxies of fitness: the number of pups produced per individual and lifespan for individuals with information on age at death. Although we lacked information about cause of deaths and reproduction is managed by humans in the ex situ population through assessments of pedigree‐relatedness using the studbook, we investigated the genomic structure of variation and autozygosity as they relate to the overall number of documented pups and age at death. We further investigated how genetic load may have responded to sharp changes in demography (i.e., inter‐lineage crosses) and whether it was associated with fitness‐related traits. We annotated the predicted protein functional impact for each locus using Ensembl Variant Effect Predictor (VEP; McLaren et al. 2016). We focused on loci that contained alleles with a predicted damage score of SIFT ≤ 0.01, which indicate a predicted amino acid substitution that would impact protein function (McLaren et al. 2016). We tabulated the number of heterozygous and homozygous genotypes per wolf for these loci of interest. We interpreted the number of homozygous loci as the realised (expressed) genetic load while the number of heterozygous genotypes as the masked (inbreeding and potential) genetic load, following theory that realised load is expected for recessive deleterious alleles expressed in homozygous genotypes while the masked load are the recessive alleles carried in heterozygous genotypes (Bertorelle et al. 2022).

Statistical Significance

2.8

Unless noted, all inter‐group differences were assessed using Welch's two‐sample t‐test, a one‐way ANOVA, Tukey HSD post hoc test, linear regression and Pearson's correlation to test for significance (p < 0.05) in R.

Results

3

We collected genome‐wide genotype data from RADseq for 187 pedigreed Mexican grey wolves (Figure 1A). Prior to data processing and filtering, we obtained 342,843 loci with an average sequencing depth of 19.5‐fold coverage (SD = 10.6, min = 5.5) from fragments with a mean insert length of 318.6 nucleotides (SD = 90.3). We retained 179 samples that did not have an excess of missing genotypes for the 134,175 loci that comprised the minimally filtered dataset, with a genotyping rate of 94.2% (8822 loci excluded due to missingness; 36,888 loci removed due to MAF threshold; eight individuals were removed due to missingness) (Table S1). Our demographic dataset, considered statistically neutral and unlinked, had 4282 loci after removing 129,893 loci that exceeded the genotypic correlation threshold and 667 loci that deviated from equilibrium (genotyping rate of 93.6%).

Breeding Population of Mexican Grey Wolves Has Declined Over the Past 50 Generations

3.1

Genetic effective population size (N e) inference from the breeding population revealed that Mexican grey wolves have significantly declined over the past 50 generations (Spearman's rank correlation: ρ = 0.991, S = 116, p < 2.20 × 10^−16^; β = 4.7152, SE = 0.35, t = 13.6, R ^2^ = 0.82, p < 2.20 × 10^−16^), from N e = 429.1 to 24.5, respectively (Figure 1B). A log transformation of the N e reveals an average decrease in N e of 2.1% per generation (β = 2.1, SE = 0.004, t = 5.60, R ^2^ = 0.43, p < 1.57 × 10^−6^). Despite protections under the ESA, this population experienced a recent decline since their listing in 1976 (pre‐1976 N e = 331.3 [324.1, 359.6], post‐1976 N e = 240.1 [226.0, 263.3]; Welch's t‐test: t = −4.22, df = 15.04, p = 0.0007). Further, we report that effective population size continued to decline after the inter‐lineage cross policy was implemented in 1995 (~21–14 generations ago) (pre‐1995 N e = 322.5 [319.6, 335.3], post‐1995 N e = 220.0 [201.1, 244.9]; Welch's t‐test: t = −3.0, df = 6.7, p = 0.0202). However, present‐day effective size estimates (mean N e between generations 1 and 13 = 211.8) are low, albeit not significantly different from the apparent stabilised N e immediately following the ESA (generations 14 and 21 N e = 269.7; Welch's t‐test: t = −1.53, df = 5.0, p = 0.1855).

Life History and Inbreeding Patterns in Mexican Grey Wolves

3.2

A PCA of neutral sites (percent variance: PC = 8.6%, PC2 = 4.1%) revealed that metrics of inbreeding significantly explained genetic variation along the first principal component (F: β = 0.583, SE = 0.16, t = 3.61, R ^2^ = 0.0685, p = 0.004; F ROH: β = 1.600, SE = 0.19, t = 8.24, R ^2^ = 0.277, p = 3.66 × 10^−14^), whereas birth year significantly described genetic variation along the second principal component (β = 0.0063, SE = 0.002, t = 4.08, R ^2^ = 0.086, p = 6.83 × 10^−5^) (Figure S1). These trends remained consistent when comparing genetic variation on an MDS analysis of identity by descent probability (IBDp) values (F on C1: β = −0.160, SE = 0.04, t = 3.71, R ^2^ = 0.072, p = 2.83 × 10^−4^; F ROH on C1: β = −0.424, SE = 0.05, t = −8.15, R ^2^ = 0.273, p < 6.37 × 10^−14^, birth year on C2: β = −0.001, SE = 0.0004, t = −3.56, R ^2^ = 0.0668, p = 4.77 × 10^−4^) (Figure S2).

Exploring patterns of inbreeding for all individuals in a given year, we observed a significant positive correlation between time and F (R = 0.835, p = 1.70 × 10^−5^), F ROH (R = 0.842, p = 1.20 × 10^−5^), number of segments in a run of homozygosity (R = 0.79, p = 9.60 × 10^−5^) and the length of runs of homozygosity (R = 0.72, p = 7.70 × 10^−5^) from the beginning of the ex situ breeding program to 1995 (Figure 1C). Following inter‐lineage crossings in 1995, we observed a dampening of inbreeding estimates through time for F (R = 0.13, p = 0.518; linear model, year × period interaction: β = 0.0146, SE = 0.00213, t = 6.88, p = 2.45 × 10^−8^), F ROH (R = −0.035, p = 0.863; linear model, year × period interaction: β = 0.0098, SE = 0.00148, t = 6.61, p = 5.88 × 10^−8^), number of segments in a run of homozygosity (R = 0.525, p = 0.005; linear model, year × period interaction: β = 0.866, SE = 0.227, t = 3.81, p = 4.53 × 10^−4^) and the length of runs of homozygosity (R = −0.266, p = 0.179; linear model, year × period interaction: β = 235, SE = 43, t = 5.46, p = 2.52 × 10^−6^) (Figure 1C). To examine relationships among life history traits, we excluded surviving individuals not yet able to reach their total reproductive output (n = 113), and found fecundity traits like number of litters and number of pups were strongly correlated with one another (R = 0.67, p = 3.10 × 10^−16^), with a significant positive correlation between lifespan and fecundity (total number of litters: R = 0.41, p = 8.2 × 10^−6^; total number of pups: R = 0.39, p = 2.0 × 10^−5^) (Figure S3). When comparing life history traits with metrics of inbreeding for these non‐surviving individuals, we found that lifespan negatively correlated with observed homozygosity (R = −0.21, p = 0.027) and the number of segments in a run of homozygosity (R = −0.2, p = 0.031) (Figure 1D). Additionally, the number of segments in a run of homozygosity (R = −0.21, p = 0.029) (Figure 1E), and the length of those segments (R = −0.27, p = 0.0043) (Figure 1F) declined with the number of pups produced by an individual.

Each Ex Situ Founder Lineages Are a Discrete Genetic Cluster

3.3

We assigned 79 query individuals to Aragón, 72 to Ghost Ranch and 28 to McBride lineages following their posterior probability population assignment values using the Bayesian clustering model STRUCTURE. The PCA clusters were concordant with founder lineage memberships, with the first principal component capturing genomic variation across all lineages (one‐way ANOVA: df = 2, F = 106.3, p < 2.0 × 10^−16^), and the second significantly distinguishing Ghost Ranch (df = 2, F = 69.6, p < 2.0 × 10^−16^) (Figure 2A). Similarly, multidimensional scaling of IBDp with respect to founder lineage membership also supported the clear and significant distinction among founder lineages on the first component of variation (one‐way ANOVA: df = 2, F = 95.3, p < 2.0 × 10^−16^), and the second distinguishing Ghost Ranch (df = 2, F = 79.4, p < 2.0 × 10^−16^) (Figure S4).

(A) PCA of unlinked SNPs, with points coloured by maximum inferred ancestry assignment to the founder lineages. (B) Inbreeding coefficient (F ROH), (C) distribution of run‐of‐homozygosity lengths (Kb) and (D) decay of linkage disequilibrium for each founder lineage assignment (standard error depicted as the shading around the mean line).

Highest Amounts of Inbreeding and Autozygosity Found in the McBride Lineage

3.4

We found a significant difference in homozygosity collectively across founder lineages, with the highest level of observed homozygosity found in the McBride lineage (mean H O, Aragón = 0.601, Ghost Ranch = 0.608, McBride = 0.637; one‐way ANOVA: F = 7.011, df = 2, p = 0.0012). Estimates of inbreeding (F) were also significantly different across founder lineages, with McBride containing the only non‐zero estimate of inbreeding (mean F, McBride = 0.023, Aragón = −0.074, Ghost Ranch = −0.055; one‐way ANOVA: F = 6.71, df = 2, p = 0.0016). Aragón and Ghost Ranch were not significantly different from each other in observed homozygosity (Tukey's HSD p = 0.6458), nor in estimates of F (Tukey's HSD p = 0.6100). However, all lineages differed in their mean estimates of F ROH, with a larger value in McBride (mean F ROH, Aragón = 0.170, Ghost Ranch = 0.207, McBride = 0.282; one‐way ANOVA: F = 18.7, df = 2, p = 4.25 × 10^−8^) (Figure 2B) and, on average, longer tracts of autozygosity than in Aragón or Ghost Ranch (mean [median] ROH size in Mb: McBride = 621 [637], Aragón = 374 [358], Ghost Ranch = 456 [509]; one‐way ANOVA: F = 18.72, df = 2, p = 4.25 × 10^−16^) (Figure 2C; Table 2). After binning inter‐SNP distances between sites in LD by 1000 Kb, we found a significant difference across all lineages (one‐way ANOVA: F = 23,661, df = 2, p < 2.0 × 10^−16^), with a higher mean r ^2^ in McBride (r ^2^ = 0.772), followed by Ghost Ranch (r ^2^ = 0.695) and Aragón (r ^2^ = 0.649) (Figure 2D).

Life History and Inbreeding Patterns Within Founder Lineages

3.5

We estimated the ancestry proportions contributed by each founder lineage prior to the inter‐lineage crossings and found significant differences in admixture across the posterior‐probability‐informed ancestry assignments (permanova df = 2, R ^2^ = 0.498, F = 87.248, p = 0.001) (Figure S5). We removed one ancestry component to account for compositional dependence and found that all groups differed significantly after FDR correction with a post hoc pairwise permutation MANOVA. This indicates that the posterior‐probability founder lineage assignments reflected distinct ancestry profiles (adjusted p = 0.001 for all comparisons; McBride vs. Aragón: F = 122.96, R ^2^ = 0.539; McBride vs. Ghost Ranch: F = 66.70, R ^2^ = 0.405; Aragón vs. Ghost Ranch: F = 74.79, R ^2^ = 0.334) (Figure 3A).

(A) Triangle plot of the supervised admixture assignment probability (Q) for each query individual with respect to the three reference founder lineages (posterior probability lineage assignment indicated by colour). (B) The number of breeding individuals assigned to each lineage ancestry and the (C) mean measures of autozygosity (F ROH) for each assigned lineage ancestry within the dataset across each year (vertical dashed line marks the 1995 inter‐lineage crossing start date). (D) Average number of heterozygous genotypes with predicted deleterious variants across each assigned lineage ancestry.

Exploring patterns for all breeding individuals in a given year, we found that the number of individuals assigned to each founder lineage increased significantly across generations for Aragón (R = 0.939, p = 1.33 × 10^−21^) and Ghost Ranch (R = 0.963, p = 2.53 × 10^−24^), whereas McBride showed no significant growth (R = 0.0725, p = 0.661) (Figure 3B). Following the inter‐lineage crossings in 1995, we found that the number of ancestry assignments through time was strongly negative for McBride (R = −0.944, p = −4.7 × 10^−14^), while strongly positive for Ghost Ranch (R = 0.783, p = 8.4 × 10^−7^) and Aragón (R = 0.840, p = 2.3 × 10^−8^) (Figure 3B). We also tested whether founder ancestry shaped patterns of autozygosity (F ROH) across time, particularly before and after the 1995 lineage intercross. A mixed‐effects model with F ROH as the response, ancestry assignment × period as fixed effects and a random intercept for year to account for variation in the number of individuals among sampling years revealed a significant interaction (F(2,77.64) = 31.40, p = 1.02 × 10^−10^). Post‐1995 F ROH differed significantly among ancestry groups (one‐way ANOVA: df = 2, F = 00.5, p = 2.0 × 10^−16^). Post hoc Tukey HSD comparisons showed higher autozygosity in McBride than Ghost Ranch (diff = 0.101, 95% CI = 0.0791–0.123, p < 0.001) and Aragón (diff = 0.121, 95% CI = 0.0991–0.143, p < 0.001), while Aragón and Ghost Ranch did not differ significantly (difference = −0.02, 95% CI = −0.042 to −0.002, p = 0.08) (Figure 3C).

After excluding surviving individuals not yet able to reach their total reproductive output (n = 113), we observed similar patterns of life history traits across Aragón and Ghost Ranch lineages, with lifespan significantly influencing the number of pups (Aragón: R = 0.317, p = 0.0264; Ghost Ranch: R = 0.551, p = 2.33 × 10^−4^) and number of litters observed throughout a lifetime (Aragón: R = 0.435, p = 0.0018; Ghost Ranch: R = 0.508, p = 8.21 × 10^−4^) (Figure S6). However, only Ghost Ranch exhibited significant negative correlations between the number of homozygous segments and lifespan (R = −0.350, p = 0.0266), as well as the number of pups (R = −0.342, p = 0.0308) (Figure S7). Additionally, the number of pups for Ghost Ranch was negatively correlated with the average length of runs in homozygosity (R = −0.402, p = 0.0101) (Figure S7).

Realised Genetic Load Was Rare but Masked Genetic Load Is Present

3.6

We tabulated the number of loci containing alleles in each of the predicted functional impact categories (Table S2). We focused on the subset of 1328 loci that have a conservative and high confidence in the putatively damaging score (SIFT ≤ 0.01), with which 263 were predicted to be highly damaging. We found no significant difference in realised genetic load between founder ancestries (mean ± SD: Aragón = 43.6 ± 17.8, Ghost Ranch = 41.2 ± 16.4, McBride = 41.9 ± 9.2; one‐way ANOVA: df = 2, F = 0.459, p = 0.633); however, Aragón carried significantly more heterozygous loci (209 ± 28.5) relative to the other lineages (Ghost Ranch = 197.1 ± 33.9, McBride = 171.1 ± 18.9; df = 2, F = 17, p = 1.78 × 10^−7^; TukeyHSD: Aragón‐Ghost Ranch p = 0.0386, Aragón‐McBride p = 1.00 × 10^−7^, Ghost Ranch‐McBride p = 3.48 × 10^−4^) (Figure 3D). We found no evidence of significant correlations between genetic load and reproduction or lifespan for each founder lineage, except for the marginally significant positive trend in Ghost Ranch between pup production and realised genetic load (R = 0.22, p = 0.058).

We found evidence that wolves born since the implementation of inter‐lineage cross policy carry significantly higher levels of masked genetic load than the founders (mean Founders = 175.8, Inter‐lineage = 202.3, 1‐tailed t‐test of unequal variance p = 1.11 × 10^−5^), which carried a higher realised genetic load albeit non‐significant (Founders = 47.9, Inter‐lineage = 41.4, p = 0.0829). We found that most assessments of genetic load over time and as a correlate of fitness were non‐significant, with only a significant positive trend found of masked genetic load and lifespan for wolves born from inter‐lineage crosses (R = 0.23, p = 0.032). One limitation is that the founder samples analysed here had an underrepresentation of Aragón and Ghost Ranch, thus results are most relevant for understanding genetic load in the McBride lineage.

Of the 263 loci predicted to have high functional impact, only three of those were associated with an annotated gene: GOSR2 (CFA9:10086828, associated with epilepsy and muscular dystrophy), KASH5 (CFA1:106678920, associated with premature ovarian and spermatogenic failure) and FEZF1 (CFA14:59519738, associated with hypogonadotropic hypogonadism). We found that homozygotes were rare with no frequency differences between founder lineages (mean ± s.d.: Aragón = 0.33 ± 0.5, Ghost Ranch = 0.28 ± 0.5, McBride = 0.25 ± 0.5; 1‐tailed t‐test of unequal variance: p > 0.05), while Ghost Ranch was significantly least likely to carry a putatively damaging allele in the heterozygous state (Aragón = 1.33 ± 0.8, Ghost Ranch = 0.96 ± 0.8, McBride = 1.32 ± 0.9; 1‐tailed t‐test of unequal variance: Aragón‐Ghost Ranch p = 0.0029, Aragón‐McBride p = 0.4848, Ghost Ranch‐McBride p = 0.0401).

Discussion

4

Following their extirpation in the wild, conservation strategies for the ex situ endangered Mexican grey wolf subspecies have relied on ex situ pedigree management driven by evolutionary principles promoting population growth and genetic diversity (Fredrickson et al. 2007; Hedrick and Fredrickson 2008; Carroll et al. 2014; Harding et al. 2016). Originally established with a single certified lineage of three unrelated individuals (McBride), there was a clear need to incorporate new genetic diversity into the severely bottlenecked pedigree (García‐Moreno et al. 1996; Harding et al. 2016). Two individuals from each of two subsequently discovered lineages (Aragón and Ghost Ranch) increased the overall genetic variation catalogued for Mexican wolves (García‐Moreno et al. 1996). However, the three lineages were managed as isolated populations until 1995, when these two additional lineages were genetically certified and approved for inter‐lineage pairings to dampen the elevated levels of inbreeding suspected from pedigree‐based analyses (García‐Moreno et al. 1996; Hedrick and Fredrickson 2008). The success of the inter‐lineage crossing immediately combated evidence of inbreeding depression, which was empirically documented through the proportion of stillbirths, litter size and pup survival (Fredrickson et al. 2007; Hedrick and Fredrickson 2008). This crossbreeding program thus began the era of genetic management for Mexican grey wolf conservation planning (Miller 2006). Although pedigree‐based population viability models provide the best available science to guide discussions surrounding genetic management, there is a pressing need for genome‐wide studies to directly evaluate the viability of the ex situ populations. As such, our genome‐wide survey of the endangered Mexican grey wolf ex situ managed population provides compelling evidence of demographic decline and marked vulnerabilities among founder lineages, offering critical insights that can assist ex situ breeding strategies in supporting genetic diversity and long‐term persistence.

Our inference of recent effective population sizes within the subset of breeding individuals (Figure 1A) revealed a steep and sustained decline over the past 50 generations (100–300 years before present), collapsing from more than 400 wolves to fewer than 25 individuals (Figure 1B). Similarly, the studbook pedigree as a whole has documented a significant decrease in the number of pups produced since its establishment in 1976 (R = −0.17, p = 2.1 × 10^−15^). However, our analyses observed a potential stabilisation of effective population size after protections under the ESA in 1976 and following the start of inter‐lineage crossing in 1995. Yet, ex situ managed breeding is constrained by a strict carrying capacity, dictated by finite space, resources and as a result of wolf releases into the wild each year that may limit or expand the number of breeding pairs possible. For example, the ex situ population experienced a slight increase in pup production in 2010, roughly 3.5 generations after the first inter‐lineage reproduction occurred, which was likely due to the cessation of moving ex situ wolves to the in situ population, allowing for more breeding pairs (Wayne and Hedrick 2010; Harding et al. 2016). Although the current breeding program continues to mitigate immediate risks of extinction, fully arresting genetic erosion is unlikely to occur.

The ex situ breeding population generally exhibited an increase in inbreeding and autozygosity from the onset of managed breeding, although this trend appears to have dampened following genetic enrichment from inter‐lineage crosses (Figure 1C). Within this population, we found a significant negative association between measures of autozygosity and fitness, with less‐inbred individuals tending to live longer and produce more offspring (Figure 1D,F). Although the ex situ population shows limited genetic structure, likely due to a small number of founders and pedigree‐informed management efforts to maintain diverse inter‐lineage pairings, these trends are likely driven by genomic variation shaped by founder ancestry.

Indeed, founder ancestry remains a strong influence on genomic variation within the population, including descendants from inter‐lineage crosses (Figure 2A). Based on estimates of ancestry, individuals assigned to the McBride founder lineage exhibited higher genomic estimates of inbreeding and autozygosity than individuals assigned to the Ghost Ranch or Aragón lineages (Figure 2B). Concomitantly, higher inbreeding was structured across the genome as longer tracts of autozygosity and statistical linkage disequilibrium, with the McBride lineage carrying the largest autozygous blocks in the highest relative abundance among the founder lineages (Figure 2C,D). In comparison, wolves carrying predominantly Ghost Ranch or Aragón founder ancestry were found to have significant fitness consequences of inbreeding on lifespan and fecundity, with higher measures of autozygosity significantly associated with shorter lifespans and lower reproductive output, patterns that may be driving trends observed across the entire sampled breeding population.

The observed population genetic patterns are likely shaped by significant founder lineage admixture in individuals from Ghost Ranch and Aragón lineages (Figure 3A). Although the McBride lineage was initially founded by three unrelated individuals, compared with the Ghost Ranch and Aragón lineages which were each founded by two, extensive inbreeding before inter‐lineage crosses began in 1995 was thought to have led to severe homogenisation of its genetic pool (García‐Moreno et al. 1996; Hedrick and Fredrickson 2008). As a result, management strategies have successfully incorporated and increased Ghost Ranch and Aragón ancestry among breeding individuals following the initiation of inter‐lineage crosses in 1995 (Figure 3B). Consistently, estimates of autozygosity (F ROH) were strongly influenced by ancestry group over time, with McBride lineages maintaining higher levels among the founder lineages following the 1995 inter‐lineage crosses despite the increased proportions of ancestry from Aragón and Ghost Ranch in the genomic background (Figure 3C).

Inbreeding depression likely manifests as reduced survival and cumulative reproductive success, underscoring concerns that genetic erosion may contribute to declines in individual fitness within the ex situ population. Our analyses revealed that patterns of genetic load differed among founder lineages. Notably, the Aragón lineage maintained significantly higher levels of masked genetic load, suggesting a potential buffering effect of standing variation in this lineage (Figure 3D). Wolves within the McBride lineage exhibited an increase in realised load, yet overall correlations with demographic performance were generally weak. However, genetic load generally decreased over time, with the inter‐lineage cross policy appearing to have successfully reduced the previously concerningly high levels of genetic load among the founders. This scenario has been previously documented in wild inbred wolf populations that experienced temporary genetic relief due to immigrants reducing the realised genetic load within the population (Smeds and Ellegren 2022). Additionally, wolves born since its implementation carried significantly higher levels of masked load compared to the founders, suggesting that recombination across divergent ancestries may have increased the pool of heterozygous deleterious alleles, some of which could remain hidden from selection. Yet, whether the concern is manifested or reflects a masked genetic load, either scenario can compromise the long‐term fitness and persistence of the ex situ population, particularly when inbreeding remains elevated following a demographic bottleneck (Dussex et al. 2023; Grossen and Ramakrishnan 2024; Jensen et al. 2024; Speak et al. 2024).

Among grey wolves, the Mexican grey wolf population exhibits a level of genetic diversity that is not as eroded as the small, inbred population of the Isle Royale National Park grey wolves following long‐term isolation (mean autosomal F ROH, ex situ Mexican = 0.20, Isle Royale = 0.35 estimated from 15 wolves; vonHoldt, DeCandia, et al. 2024, vonHoldt, Stahler, et al. 2024). However, Mexican grey wolves exhibit lower autozygosity in comparison to reintroduced grey wolves of the Northern Rocky Mountains (F ROH = 0.31 estimated from 136 wolves), as well as those in the Pacific Northwest (F ROH = 0.27 estimated from 38 wolves), experiencing recent founder effects from similar gene pools and small effective population sizes maintained at roughly 9% of the census estimates (vonHoldt, DeCandia, et al. 2024, vonHoldt, Stahler, et al. 2024). Similarly, the small isolated grey wolf population of the Sierra Morena on the Iberian Peninsula has also persisted via few packs with high levels of autozygosity, with genomic analyses revealing recent introgression from domestic dogs drifting to increased frequencies within the population, threatening large genomic fragments of wolf ancestry with replacement and ultimate extinction (Gómez‐Sánchez et al. 2018).

To predict the long‐term conservation implications of low genetic diversity in the severely bottlenecked Mexican grey wolf population, we can draw on patterns of inbreeding and genetic load documented in other systems, particularly endangered species that have undergone genetic rescue and restoration (Hedrick and Fredrickson 2010). One analogous system to the Mexican grey wolf is the history and conservation management of the Pacific pocket mouse ( Perognathus longimembris pacificus ). Historically found along southern California's coast, this endemic species was thought to be extinct by 1935 (Patton and Fischer 2023; Williams 1986). In 1993, a single population was rediscovered and listed as endangered (Dana Point Headlands), with two more populations identified (North and South San Mateo) in 1995 (USFWS 1998; Brylski et al. 1998; Spencer et al. 2000; Spencer 2005). Yet wild populations continue to face challenges as they are surrounded by an urban landscape, thwarting dispersal and gene flow, with < 50 individuals reportedly surviving at each site (Brylski et al. 1998; Swei et al. 2003; Spencer 2005). Consequently, an ex situ breeding program was established at the San Diego Zoo Institute for Conservation Research in 2012 (Shier et al. 2016; Chock et al. 2018). Although this population was founded by 10 individuals from each of the three extant lineages, it was managed as a single population due to their limited evolutionary divergence (Shier et al. 2016; Chock et al. 2018). This inter‐lineage reproduction initially increased fitness, consistent with genetic rescue, yet evidence has suggested that this early benefit may mask potential long‐term risks (Wilder et al. 2022). Recent genomic assessments suggest possible elevations in genetic load associated with founder‐lineage ancestry (King et al. 2022), resulting from inbreeding and low effective population sizes (Wilder et al. 2020, 2022; Lynch 1991; Hedrick and García‐Dorado 2016). Effective management therefore requires strategies that minimise realised load through the reduction of inbreeding, while also curbing masked load by enforcing stringent selection practices that ensure only individuals carrying the lowest levels of deleterious variation are prioritised for breeding or release.

Such conservation efforts have long worked to restore gene flow for demographic and genetic rescue of other endangered species. Some of the most classic examples include the endangered subspecies Florida panther ( Puma concolor coryi ) genetic rescue program that was initiated in 1995 for intentional targeted introgression of 25% ancestry from their proximal Texas subspecies (Pimm et al. 2006; Johnson et al. 2010; Hostetler et al. 2013; Van De Kerk et al. 2019; Onorato et al. 2024; Aguilar‐Gómez et al. 2025) or the greater prairie chicken ( Tympanuchus cupido pinnatus ) demographic rescue via translocations between 1992 and 1998 (Svedarsky et al. 2000; Mussmann et al. 2017). Classic examples of documented genetic rescue via natural dispersal that have demonstrated how new genetic variation can relieve the detrimental effects of deleterious genetic load include Isle Royale (Wayne et al. 1991; Adams et al. 2011; Hedrick et al. 2014; Robinson et al. 2019) and Scandinavian grey wolves (Vilà et al. 2003; Kardos et al. 2018).

While the ESA was enacted in 1973, Mexican grey wolves were not listed for protection until 1976. At that time, wolves were documented to be extirpated from the wild in the United States, and only 80 wolves were estimated to remain in Mexico (per. comm., R. McBride 1980). Although this population decline is known to have resulted in low levels of genomic diversity and high inbreeding within remnant populations (Leonard et al. 2005; Taron et al. 2021), questions remain regarding best management practices for the Mexican grey wolf subspecies. Even though the U.S. and Mexico ex‐situ populations are technically managed as a single demographic unit, logistical challenges can limit cross‐border transfers. As a result, Mexico is suspected to have distinct differences in lineage ancestry due to different compositions of founder individuals and breeding strategies. Therefore, broader international population genomic surveillance is needed to further explore how future management decisions could support increased retention of genetic diversity within the subspecies, particularly maintaining valuable genetic contributions from the Ghost Ranch and Aragón lineages. Genomic restoration of historical diversity would be significantly difficult; however, potential recovery of genomic diversity from historical museum specimens may eventually help revive lost variation in this subspecies via unconventional biotechnological advances (Fritts 2022; Howard et al. 2016). Since the Mexican grey wolf recovery plan was initiated, federal support has worked to carefully archive all ex situ and in situ wolves in a permanent repository to provide the biomaterials necessary for such revolutionary efforts, as well as stimulate investigations into other aspects of wolf biology related to their conservation (e.g., National Academy of Sciences, Engineering, and Medicine 2021; Witt et al. 2024). Ultimately, successful conservation will require more than restoring lost genomic variation; it will depend on the collective willingness to confront unresolved management divides, secure sustained resources and meaningfully engage the communities most affected by wolf recovery. This specifically rests upon refreshing a coordinated binational management and integrating Indigenous knowledge with local community partnerships. All are essential for securing the resources and public support needed to restore a genetically resilient and socially supported wild Mexican grey wolf population.

Author Contributions

B.M.v. and Y.L. designed the research; M.L.C., J.L.D. and J.A.C. prepared, archived and digitised specimens over a 20 year period; Y.L., M.K., E.K., I.G.N., Z.W. and B.M.v. performed the research analysis; all authors contributed to writing and revising the manuscript.

Funding

This study was supported by Princeton University's F250th educational development award to B.M.v.

Disclosure

Benefits Sharing: Benefits from this research accrue from the sharing of our data and results with federal and state agencies and public repositories. More broadly, our research explores questions about the life histories and demographic patterns that likely impacts the conservation of an endangered species. We are committed to building institutional infrastructure and human capacity that will enable the preservation of genetic diversity and the future persistence of an endangered species, with a focus on binational relationship building and implementation of recovery plans.

Conflicts of Interest

B.M. is an active member of the Red Wolf Recovery Team and has received funding from the US Fish and Wildlife Service for genomic studies of red wolves and coyotes along the Gulf Coast of the United States but declares no conflicts of interest with respect to this paper.

Supporting information

Figure S1: Determinants of variation within a PCA from unlinked SNPs, including sex, year of birth, age at death, number of litters (top, left to right), pedigree inbreeding estimates (F), inbreeding estimates from autozygosity (F ROH), number of segments in a run of homozygosity and the average length (Kb) of a segment in a run of homozygosity (bottom, left to right). Figure S2: Determinants of variation within an MDS plot of identity by descent probability (IBDp) from unlinked SNPs, including sex, year of birth, age at death, number of litters (top, left to right), pedigree inbreeding estimates (F), inbreeding estimates from autozygosity (F ROH), number of segments in a run of homozygosity and the average length (Kb) of a segment in a run of homozygosity (bottom, left to right). Figure S3: Relationships between fecundity (left) and lifespan (centre, right) across all non‐surviving individuals in the ex situ population genomic dataset. Linear model with 95% confidence intervals shown for all panels. Figure S4: Determinants of variation within an MDS plot of identity by descent probability (IBDp) from unlinked SNPs, with points coloured by maximum inferred ancestry assignment to the founder lineages. Figure S5: Supervised admixture ancestry probability (Q) for each Mexican grey wolf grouped by their maximum posterior probability lineage assignments. Figure S6: Relationships between life history traits within assigned founder lineage assignment (loess‐smoothed linear model with 95% confidence intervals provided). Figure S7: Relationship between the number of pups and run of homozygosity (ROH) segment length (left), the number of runs of homozygosity (centre) and the number of runs of homozygosity and lifespan (right) within the Ghost Ranch lineage. Linear model with 95% confidence intervals shown for all panels. Table S2: Annotation of loci out of a total of 267,892 loci.

Table S1: Mexican grey wolf meta‐data as it pertains to studbook life history details (sex, pedigree parents, lifespan, birth location, number of pups). Individuals included in this analysis are indicated by the CL# as the identification number of the DNA sample (F, inbreeding coefficient; F ROH, inbreeding coefficient derived from runs of homozygosity across the autosomes; H E, expected heterozygosity; H O, observed heterozygosity; K, genetic cluster; Q, posterior genetic cluster assignment probability; YOB, year of birth; YOD, year of death).