Genomic recovery lags behind demographic recovery in bottlenecked populations of the Channel Island fox, Urocyon littoralis

With continued global change, recovery of species listed under the Endangered Species Act is increasingly challenging. One rare success was the recovery and delisting of the Channel Island fox (Urocyon littoralis) after 90%–99% population declines in the 1990s. While their demographic recovery was marked, less is known about their genetic recovery. To address genetic changes, we conducted the first multi‐individual and population‐level direct genetic comparison of samples collected before and after the recent bottlenecks. Using whole‐exome sequencing, we found that already genetically depauperate populations were further degraded by the 1990s declines and remain low, particularly on San Miguel and Santa Rosa Islands, which underwent the most severe bottlenecks. The two other islands that experienced recent bottlenecks (Santa Cruz, and Santa Catalina islands) showed mixed results based on multiple metrics of genetic diversity. Previous island fox genomics studies showed low genetic diversity before the declines and no change after the demographic recovery, thus this is the first study to show a decrease in genetic diversity over time in U. littoralis. Additionally, we found that divergence between populations consistently increased over time, complicating prospects for using inter‐island translocation as a conservation tool. The Santa Catalina subspecies is now federally listed as threatened, yet other de‐listed subspecies are still recovering genetic variation which may limit their ability to adapt to changing environmental conditions. This study further demonstrates that species conservation is more complex than population size and that some island fox populations are not yet ‘out of the woods’.


| INTRODUC TI ON
The recovery of endangered species is a rare occurrence. Since passage of the Endangered Species Act (ESA) in 1973, only 39 out of 2300 species were delisted due to recovery as of 2017 (U. S. Fish and Wildlife Service, 2017). A report on endangered birds estimated that recovery plans would take an average of 63 years (Suckling et al., 2016), illustrating that recovery is generally a slow process. In contrast, the Channel Island fox, Urocyon littoralis, was delisted after spending just 12 years on the federal endangered species list securing the record of the fastest successful recovery for any ESA-listed mammal in the United States (U.S. Department of the Interior, 2016).
Four of the six island fox subspecies were federally listed in 2004 following 90%-99% population crashes that occurred in the late 1990s on Santa Catalina and the northern Channel Islands due to canine distemper virus and hyperpredation by golden eagles respectively (Coonan et al., 2005(Coonan et al., , 2010(Coonan et al., , 2014King et al., 2014;Roemer et al., 2001;Timm et al., 2009). The two other fox subspecies, on San Clemente and San Nicolas, are not known to have undergone recent bottlenecks. Recovery on the affected islands was made possible by a combination of canine distemper vaccination, golden eagle relocation, removal of feral pigs (prey for golden eagles), re-introduction of bald eagles (competitors for golden eagles that are less likely to prey on foxes) and highly effective captive breeding programmes (Coonan et al., 2014;King et al., 2014;Roemer et al., 2001). With these efforts, the affected fox populations rebounded to near, at, or above pre-crash population census sizes (Table 1) and therefore were delisted in 2016, except for the foxes on Santa Catalina which were downlisted to threatened. On San Miguel and Santa Rosa islands, all foxes were taken into captivity making these subspecies extinct in the wild at the beginning of the captive breeding programme (Coonan et al., 2010). Meanwhile, on Santa Catalina and Santa Cruz islands, only subsets of the remaining foxes were taken into the captive breeding programmes (Coonan et al., 2010).
While the Channel Island fox demographic recovery has been marked, much less is known about their genetic recovery. Generally, genetic diversity declines alongside population decreases, but often takes longer to recover than does census population size (Allendorf, 1986;Nei et al., 1975). Therefore, we rigorously tested, for the first time, if the island foxes recovered their genetic diversity. Previously, whole-genome sequencing from one individual per island (two from San Nicolas Island) collected before the bottlenecks showed exceedingly low genetic diversity within islands, especially for the two individuals from San Nicolas (Robinson et al., 2016).
Mitogenomes sequenced from samples collected after the bottlenecks suggested that some haplotype diversity was retained despite population declines (Hofman et al., 2015). Additionally, a RADseq study from samples collected after the bottlenecks (Funk et al., 2016) showed low nucleotide diversity but higher observed and expected heterozygosity within islands than that found by Robinson et al. (2016). A subsequent whole-genome study by Robinson et al. (2018) found no statistically significant change in genome-wide heterozygosity between samples collected in 1988 TA B L E 1 Summary of population declines and results from this study. p-values are based on Wilcoxon rank sum tests between historical and modern samples within an island. Islands that went through a recent bottleneck are indicated with an asterisk. Population minimum estimates taken from (Coonan et al., 2010(Coonan et al., , 2014King et al., 2014 (Laughrin, 1980), they found consistently low genetic diversity in samples from 1929, 1988 and 2000 (1-2 individuals per year). While previous work did not detect significant temporal genetic changes, these studies were constrained by minimal sampling pre-and post-bottleneck.
In this study, we evaluate genetic diversity across multiple individuals collected pre-and post-bottleneck in each of the six U.
littoralis subspecies. Three aspects of our study that will help us accurately test if genetic diversity has recovered are the use of exome sequencing, multiple diversity metrics and multi-individual sampling for all islands pre-and post-bottleneck. An exome capture approach was chosen to reduce sequencing human and/or microbial contaminants that can be associated with museum specimens (Bi et al., 2013;Carpenter et al., 2013) and to focus on genic regions and putatively adaptive loci. While exome sequence data does not have the benefit of revealing neutral population structure, genic regions are advantageous in that they have been shown to be better at determining population differentiation than intergenic markers (Gelfman et al., 2012;Hand et al., 2016;Matala et al., 2014;Zhan et al., 2015).
Given that popular genetic metrics assess different aspects of diversity and therefore respond differently to demographic changes, we employed multiple metrics to better diagnose if U. littoralis has fully recovered. Here we measure heterozygosity, which is expected to vary with size of the bottleneck and the rate of population growth after the contraction (Nei et al., 1975). We also measure nucleotide diversity which is expected to be more sensitive than average heterozygosity to bottlenecks due to greater reliance on rare alleles (Luikart et al., 1998;Nei et al., 1975). In addition, we measure the number and length of runs of homozygosity (ROH) which are expected to increase with bottlenecks and decrease with time since bottleneck (Ceballos, Joshi, et al., 2018). We also examine divergence, F ST , among islands over time, which is expected to increase after a bottleneck due to drift. Finally, using multiple samples from pre-and post-bottlenecks across all inhabited islands we will be able to calculate population-level statistics and assess variation among individuals.
By focusing on functional regions of the genome we hope to shed light on the genetic health of island fox populations. Previous work has shown little evidence for inbreeding depression, despite extraordinarily low genomic variation. Coonan et al. (2010) remarked at the absence of detectable inbreeding symptoms in captive populations, except perhaps for some elevated perinatal mortality. Robinson et al. (2018) found that skeletal defects linked to inbreeding in wolves and cheetahs were rare in island foxes. It is important to note, however, that a lack of evidence for inbreeding depression is not evidence for an absence of inbreeding depression (Ralls et al., 2020). We take an alternative approach to assessing genomic health, by focusing on coding regions more closely tied to fitness.
We also use this quasi-natural experiment to test if the different bottlenecks left distinct genomic signatures due to the differences in cause, severity and recovery. For example, we might expect genes related to body size, muscle movement or sight to be selected for on the northern islands where the cause of the decline was predation.
However, on Santa Catalina Island, we might expect genes related to the immune system to be under selection because the population decline was caused by a disease.
Understanding the genetic response across the subspecies of U. littoralis is necessary to assess each subspecies' capacity for resilience. High-standing genetic variation is advantageous because selection among existing mutations can happen more quickly than waiting for new mutations (Barrett & Schluter, 2008 were provided by managers or contractors of the islands (see Acknowledgements). Whole-genomic DNA was extracted using a modified protocol (Wong et al., 2007), and the end of the protocol was made to match that of the historical sample extraction.
Briefly, red blood cell (RBC) lysis buffer was added to approximately 200 μL of a clotted blood sample. Two bleached, autoclaved, flamesterilized steel beads (3.175 mm) were then added to break open cells when vortexed. Buffer, Proteinase K and 10% SDS were added to the disrupted pellet which was then incubated at 56°C overnight.
The remaining steps followed the historical sample extraction, but DNA was eluted in 100 μL molecular grade water.
Genomic DNA was quantified using an Invitrogen Qubit 2.0 Fluorometer. DNA fragment size was evaluated with an Agilent Bioanalyzer 2100 (Agilent Technologies), and samples with many large fragments (>500 bp) were sheared using a Covaris S2/E210 Focused-ultrasonicator. Large DNA fragments were sheared to a target size of 200 bp with the following parameters: 10% duty cycle, intensity level five, 200 cycle bursts, frequency sweeping mode, time 60s and a total of three cycles.

| Exome capture and sequencing
We used commercially made exome capture baits (NimbleGen's SeqCap EZ Developer) based on a previous design (Broeckx et al., 2014) for the genome of the domestic dog, Canis familiaris.

| Within-population statistics
We split the filtered VCF by island and time period to calculate within-population statistics. We calculated nucleotide diversity (π) in 10 kb non-overlapping windows with Pixy v. 1.2.7.beta1 (Korunes & Samuk, 2021;10.5281/zenodo.6551490). Since Pixy highlights the importance of incorporating missing data, including invariant sites, in calculations of π, we re-called SNPs with Freebayes --report-monomorphic per scaffold then merged the VCFs using BCFtools concat. Variants were filtered to remove indels, depth less than 5, mean depth greater than 500 and more than 25% missing data using VCFtools (--remove-indels --minDP 5 --max-meanDP 500 --max-missing 0.75). Then we created separate VCFs for invariant and variant sites using VCFtools --max-maf 0 and --mac 1, respectively, following the Pixy tutorial. Individuals with greater than 75% missing data, sites with excess observed heterozygosity (greater than 0.6), variants in linkage disequilibrium (LD) were removed from the variant sites VCF following the previously stated methods. We merged the invariant and filtered variant datasets after removing individuals from the invariant data that were not in the filtered variant data. Wilcoxon rank sum tests with continuity corrections were used to determine differences between time periods within an island. Wilcoxon effect size (r) was calculated using the rstatix package (Kassambara, 2023).

| Population differentiation
Population differentiation analyses were conducted on the quality filtered, maf filtered, LD pruned dataset. We conducted a multidimensional scaling analysis (MDS; --mds-plot -cluster) in Plink and visualized it in RStudio with R. Weighted global population differentiation (Weir & Cockerham, 1984), F ST , was calculated using VCFtools (--weir-fst-pop) pairwise by island within a time period. To further evaluate relatedness among samples, we calculated identityby-descent (IBD) in Plink (--genome). We created a phylogenetic tree using SNPhylo v. 20160204 (Lee et al., 2014) which quality filters variants, generates representative sequences from SNPs and performs multiple alignment using MUSCLE v. 3.8.31 (Edgar, 2004).
With the SNPhylo output, we calculated a pairwise distance matrix and then conducted clustering with the upgma function in the R package Phangorn (Schliep, 2011). Maximum likelihood estimates were calculated for the tree with the pml function, which was followed by optim.pml to optimize the model parameters, both from Phangorn. Finally, we performed a non-parametric bootstrap analysis with 100 replicates using bootstrap.pml. The phylogenetic tree was plotted with the R package ggtree (Yu et al., 2017).

| Bottleneck and island-associated variants
Association analyses were carried out on the quality filtered data- We visualized the quantile-quantile (qq) plots and estimated the overdispersion (λ) of p-values for each test using gcontrol2 (Devlin & Roeder, 1999) in the R package Gap v. 1.1.22 (Zhao, 2007). We then looked up the significant variants using their chromosome and position in the CanFam3 reference genome on the UCSC Genome Browser (Kent et al., 2002) to find potentially corresponding genes.

| Exome capture and sequencing
We and San Nicolas (SNI, N = 17) islands (Table S1; Figure 1). We also sequenced three technical replicates-one from SRI and two from CAT-as part of the 95 samples (Table S1), but they were not used in analyses. For historical samples, the raw total PE counts ranged from 3 million to 162 million reads with an average of over 73 million total PE counts. For modern samples, we obtained raw total PE counts from over 29 million to over 138 million with an average of over 74 million total PE counts. Historical and modern sequences mapped at an average of 61.6% and 68.7%, respectively, to the dog genome (  Figure S1). Therefore, those samples were excluded from further analyses.

| Variant calling and filtering
There was a total of 27,003,261 variants in the unfiltered dataset. Initial quality filtering resulted in 6,073,646 variants and removing sites and individuals with greater than 75% missing data resulted in 3,721,301 variants (Table S3). Of these variants, 277,360 (6.3%) were identified as synonymous by snpEff. Estimates of LD ranged from a mean R 2 of 0.28 ± 0.35 in modern SNI to 0.76 ± 0.36 in modern SCZ samples ( Figure S2). Pruning for sites in LD resulted in 318,406 variants. The dataset for population differentiation was based on 104,173 variants.

| Within-population statistics
Populations known to have experienced recent bottlenecks (SMI, SRI, SCZ, CAT) generally showed signs of decreased genetic diversity. All six islands significantly decreased in mean nucleotide diversity from historical to modern time periods (Table 1; Figure 2a), however, all Wilcoxon effect sizes ranged from 0.02 (SNI) to 0.2 (SCL), which were all categorized as small in magnitude. (r SMI = .07, r SRI = .02, r SCZ = .16, r CAT = .10, r SCL = .2, r SNI = .03).
Half of the bottlenecked islands, SMI and SRI, showed a significant decrease in individual observed heterozygosity (Figure 2b) from historical to modern time periods (p SMI = .001, p SRI = .003). The Wilcoxon effect sizes (r) for both were large (r SMI = .71, r SRI = .77).
For the remaining two bottlenecked islands, SCZ showed a nonsignificant increase, while CAT showed a non-significant decrease.
For the two islands not known to have undergone a recent bottleneck, a significant decline in observed heterozygosity between the two time periods was detected in SNI (p = .004, r = .72), while a non-significant increase was detected in SCL (p = .15). Average heterozygosity calculated in 100 Kb windows ( Figure S3) shows that heterozygosity patterns were generally consistent across the genome ( Figure S3). We conducted linear models with observed heterozygosity and year collected for each island (Table S4, Figure S4) and found significant year terms for three of the four bottlenecked islands (p SMI = 2.5 × 10 −10 , p SRI = 1.5 × 10 −3 , p SCZ = .71, p CAT = .01), and half of the non-bottlenecked islands (p SCL = .42, p SNI = 2.1 × 10 −5 ). This analysis by year ( Figure S4) demonstrates that a decline in heterozygosity began before the 1990s bottlenecks in SMI, CAT and SNI. A similar analysis could not be conducted for nucleotide diversity because it was estimated from individuals collected over multiple years.
A Wilcoxon rank sum test showed a significant increase in the number of ROH ( Figure S5a, Table 1) between historical and modern samples in SMI (p = .0008), and a marginally non-significant increase in SRI (p = .063). For the length of ROH ( Figure S5b, Table 1), there was a significant increase between historical and modern on SMI (p = .0008), SRI (p = .011) and SCL (p = .003). There was no significant difference in the average length of ROH on SCZ, CAT or SNI. We did find a significant positive relationship between the number of ROH and length of ROH on all islands ( Figure S5c).
Estimates of average Tajima's D were found to be negative for both modern and historic samples across all islands (Table S5,   Table 1). We found finite estimates of the effective population sizes (Table S6) only for CAT, where the historical estimate was 8.9 (95% CI: 1.8-infinite) and the modern estimate was 60.8 (95% CI: 9.6-infinite).

| Population differentiation
Bottlenecked islands showed temporal differentiation, but greater differentiation was found among islands. The multidimensional  Table S7). The IBD analysis of pairwise samples were all either 'other related' or unrelated.
There were no pairs of samples deemed full siblings, half siblings, or parent-offspring, although within island pairs were more related than inter-island pairs ( Figure S6). The maximum likelihood tree continued to show separation by island and time period (Figure 4).
Samples from the same island generally clustered together, although support for the CAT and SNI clades was weak, and one SMI sample clustered with SRI. The southern and northern samples formed wellsupported, distinct groups.

| Bottleneck and expansion-associated variants
We conducted association tests separately for the northern islands and CAT since the cause of the decline was the same across the three northern islands but different from the cause of the decline on CAT. The northern Channel Islands showed overdispersion in significance values from the association tests, λ = 2.2 ( Figure S7a). Islands that went through a recent bottleneck are indicated with an asterisk.
We found 117 variants after a Benjamini & Hochberg false discovery rate (FDR) correction that significantly changed between historical and modern samples. Four variants were still significant after a strict Bonferroni correction ( Figure S7b, Table 2). Three of these variants were in close proximity on chromosome 2 and were associated with the von Willebrand factor A domain containing the 5B1 (VWA5B1) gene, which is linked to a plasma glycoprotein involved in platelet adhesion (Perkins et al., 1994). The remaining variant was associated with the EIPR1 gene, which is involved in protein transport (Topalidou et al., 2016). The change in allele frequencies for each significant SNP can be found in Table 2. We did not find any SNPs that significantly changed between historical and modern samples on CAT after an FDR correction.

| DISCUSS ION
Our population-wide U. littoralis exome study showed that genetic diversity decreased, in most bottlenecked populations, in contrast to the only other study that compared pre-and post-bottleneck samples, which did not find differences (Robinson et al., 2018). Generally, our MDS and phylogeny showed that the fox samples changed in genetic composition over time. Notably, SMI, which experienced one of the most severe bottlenecks and remains relatively small (Table 1), consistently lost genetic diversity as shown by decreased heterozygosity and nucleotide diversity, as well as increased frequency and length of runs of homozygosity. In addition, SMI shows particularly marked increases in divergence from the other five islands over time.
SRI, which underwent a similarly severe bottleneck but rebounded to a larger population size, also showed decreased heterozygosity and nucleotide diversity, and increased length of ROH, but showed no change in the number of ROH and more moderate increases in divergence from other islands.
Only eight individuals on SMI and 12 on SRI were genetic founders of the respective captive breeding programmes (Coonan et al., 2010). This small number of individuals can lead to founder effects in captive populations (Holgate, 1966;Nei et al., 1975). Previous estimates for a captive breeding population to maintain genetic diversity are 20-30 individuals (Frankham et al., 2010). Additionally, captive population sizes are limited by resources-captive fox populations were kept at less than 40 individuals (Coonan et al., 2010)which keeps effective population sizes small. Strong founding effects, small effective population sizes and in this case a need to grow the population quickly often lead to increased relatedness and lower diversity, which can increase divergence between captive and wild populations (Woodworth et al., 2002). The decline in genetic diversity and increased inter-island divergence we show, especially for these two island populations, is potentially a result of both the initial bottlenecks and the captive breeding programmes, despite managers basing breeding decisions on the known genetic history (Coonan et al., 2010). In fact, our results agree with the genetic data that indicated SMI had the lowest genetic diversity within the captive populations compared to the other two northern subspecies (Coonan et al., 2010).
The two islands with less severe recent population crashes, SCZ and CAT, had inconsistent signals of decreased genetic diversity for the metrics tested. The divergence between islands increased over time in SCZ foxes and they showed a decrease in nucleotide diversity, which was consistent with our predictions, and yet they had no change in heterozygosity or the number or length of ROH. CAT foxes tended to decrease in diversity over time, as indicated by a decrease in nucleotide diversity, and they also showed an increase in divergence from other islands. However, CAT foxes neither showed a significant decrease in heterozygosity nor did the number or length of ROH change over time. SCZ and CAT populations had more individuals remaining after the crash than the other two northern F I G U R E 3 The change in Weighted Weir and Cockerham F-statistic estimates (modern F ST − historical F ST ). Darker red indicates larger divergence or populations growing further apart over the time period. Islands that went through a recent bottleneck are indicated with an asterisk.
islands, and the foxes were not extinct in the wild during the captive breeding programme like the other two northern populations, which may contribute to this difference. Moreover, these results are consistent with Nei et al. (1975) indicating that heterozygosity is less sensitive to bottlenecks than allelic diversity (in this study nucleotide diversity), which also underscores the severity of the crashes on SMI and SRI. It has also been suggested that only 10 founding individuals are needed to maintain 95% of heterozygosity in captive populations (Frankham et al., 2010), which could also contribute to the maintenance of heterozygosity we see on SCZ and CAT. Since there were more individuals left on SCZ and CAT when the captive breeding programmes were started, managers were able to choose individuals from the wild and did not have to take all foxes into captivity. When the captive-born foxes were released into the wild they could mate with wild individuals which would likely maintain higher genetic diversity than the islands where this was not an option. Our finding of a decline in genetic diversity over time in U. littoralis, which was not detected in previous work (Robinson et al., 2018), may be due to our increased sample size. We found individual variation of observed heterozygosity within an island for each time slice indicating information is lost if only one individual is sampled.
Consistent with previous studies (Funk et al., 2016;Robinson et al., 2016) (Robinson et al., 2018). Heterozygosity on CAT appeared to be declining before the recent bottleneck, and foxes on SNI, a population that is not known to have undergone a recent bottleneck, also showed a multi-decadal decline in heterozygosity. Our negative Tajima's D estimates across islands and time points were consistent with the declines in heterozygosity, which could suggest purifying selection (or population growth) in these populations prior to the 1990s. While individuals on SNI had low genetic diversity for some of the metrics and showed little variation in our tree, we found that they had intermediate levels of heterozygosity compared to other islands. Our results suggest that functional genetic variation has not flatlined in SNI foxes, which may contribute to the population's apparent avoidance of documented inbreeding depression (Coonan et al., 2010;Robinson et al., 2018) and continued success. However, there could be cryptic inbreeding depression in traits not yet examined such as total lifetime fitness (Ralls et al., 2020). It is important to note, however, that in highly inbred populations like island foxes, low variation in inbreeding among individuals may make inbreeding depression difficult to detect, even with good fitness data. Recent studies on the U. littoralis microbiome suggest that it could be another trait to evaluate for inbreeding depression, and the diversity found in the U.
littoralis microbiome could aid fox populations in responding to the changes in the environment despite low genetic diversity (Adams et al., 2021;DeCandia et al., 2020). Future work would benefit from demographic modelling on the historical and modern samples from this study, especially with SNI samples to investigate potential bottlenecks around the 1970s and forward-in-time simulations to assess if diversity observed in this study would be predicted under simulated bottlenecks.
Our maximum likelihood tree shows samples from the same island clustering together, with a distinct split between northern and southern islands, consistent with previous studies (Robinson et al., 2016(Robinson et al., , 2018. There was one instance of a historic SMI individual grouping with the SRI clade, which could indicate the possible movement of animals between islands. The movement of foxes between islands is a rare event, but there is anecdotal evidence that it has happened in the modern era. Alternative explanations include mislabeling or sample degradation (Axelsson et al., 2008), although we find little evidence for post-mortem DNA damage. Despite this possible evidence of inter-island migration in the past, our F ST results clearly show that island populations are becoming more genetically divergent over time. This increase in population divergence is relevant to the possibility of using inter-island translocation for conservation purposes (e.g. Funk et al., 2016). While translocation is often believed to be advantageous (e.g. Frankham, 2016), there is a trade-off between introducing new beneficial vs. detrimental variation Ralls et al., 2020;Von Seth et al., 2021).
The likelihood of successful genetic rescue is higher if island divergence is driven by the accumulation of deleterious mutations due to genetic drift, and lower if island divergence involves purging of deleterious mutations and/or adaptation to local conditions. The possibility of local adaptation is supported by genomic data (Funk et al., 2016), and may be strengthened by the occurrence of divergent microbiomes among island populations (Adams et al., 2021).
Overdispersion of significant values and differential variants identified on the bottlenecked islands suggest a systematic difference across the genome between time periods and may point to specific underlying physiological responses induced by the population crashes. We found highly significant variants between historical and modern samples on the northern islands for two genes: EIPR1 and VWA5B1. The EIPR1 gene is involved in protein transport in endosomes and the Golgi apparatus (Topalidou et al., 2016). Specifically, the EIPR1 gene is part of the endosome-associated recycling protein (EARP) complex and is important in the endocrine and nervous system because it is necessary for the regulation of dense-core vesicles, which contain molecules such as peptide hormones and neurotransmitters (Topalidou et al., 2016). The VWA5B1 gene has been shown to be expressed in mice in male testes and male and female pituitary glands (Naqvi et al., 2019), suggesting it plays a role in the reproductive and nervous systems. While we did not find evidence of selection in the genes that we had originally predicted to be important (those related to body size, muscle movement and sight on the northern islands and the immune system on Santa Catalina Island), we did identify novel candidates that appear to respond differently to population declines driven by different mechanisms.
The candidate genes putatively under selection could also be driven by the captive breeding programmes, which can be intense selection events Woodworth et al., 2002).
This study demonstrates that a multi-individual population-level study is required to gain nuance in the understanding of genetic consequences following severe population bottlenecks in U. littoralis. In contrast to Robinson et al. (2018), we found significant differences between time periods within islands, with particularly large drops in genetic diversity on the two islands that experienced the largest population declines. These drops are notably large, given that declines in heterozygosity are expected to be small if the population recovers from a bottleneck quickly or if the bottleneck comprises more than two individuals (Nei et al., 1975), both of which appear to be true in island foxes. Theoretically, in the absence of migration, genetic diversity should increase in a population at a rate inverse to that of the mutation rate (Nei et al., 1975). Using the dog mutation rate of 10 −8 per generation (Freedman et al., 2014;Lindblad-Toh et al., 2005) and an assumed generation time of 1 year (Coonan et al., 2010), it would take ~10 8 generations, or 100 million years for the island foxes to recover their genetic diversity. This time frame suggests that we would not expect the genetic diversity of the populations that experienced a decline to have recovered, however, all island populations were delisted, except for CAT. Our study suggests that the genetic diversity in CAT foxes may be lower than previously thought, which in addition to persistent anthropogenic threats supports its continued listing under the Endangered Species Act. Further, the dramatic demographic recovery on CAT has produced high fox densities that are likely to increase disease transmission (Sanchez & Hudgens, 2015).
Broadly, our study supports the notion that population genetics metrics respond differently to changes in demography and that using multiple metrics to assess genetic recovery is advisable. In our study, heterozygosity provided a more conservative metric of the genetic consequences of bottlenecks than did nucleotide diversity.
Moreover, the northern Channel Islands were delisted from the federal endangered species list because their census population sizes rebounded to pre-crash numbers, yet our results showed a lag in genetic diversity recovery, especially in SMI, suggesting that the island foxes may still struggle with long-term adaptation. These results join previous studies demonstrating the importance of assessing population viability using census and genetic data. For example, an extreme bottleneck in the little spotted kiwi (Apteryx owenii) decreased genetic diversity that still has not recovered despite also being considered a 'conservation success' and being downgraded in conservation status (Taylor et al., 2017). Similar stories of downlisted species that were demographically recovered but still had low genetic diversity are found in the dawn redwood tree (Metasequoia glyptostroboides; Li et al., 2012), the Seychelles Magpie-robin (Copsychus sechellarum; Cavill et al., 2022) and the pink pigeon (Nesoenas mayeri ;Jackson et al., 2022). Better incorporation of genetic data, including longterm monitoring, in threatened and endangered species listings will improve management of these populations that are less able to respond to novel challenges.

AUTH O R S ' CO NTR I B UTI O N S
Suzanne Edmands and Nicole E. Adams designed the study. Nicole E.
Adams performed the laboratory work and analysed the data. Nicole E. Adams and Suzanne Edmands wrote the manuscript.

ACK N O WLE D G E M ENTS
We thank the following people and institutions for their sample con-

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors have no conflict of interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
Raw exome sequence reads and metadata are available on NCBI's Sequence Read Archive digital repository (BioProject PRJNA970412), and R scripts are available on GitHub (https:// github.com/Nicol eAdam s-sci/islan dFox_exome).