Genomics Reveals Complex Population History and Unexpected Diversity of Eurasian Otters (Lutra lutra) in Britain Relative to Genetic Methods

Abstract Conservation genetic analyses of many endangered species have been based on genotyping of microsatellite loci and sequencing of short fragments of mtDNA. The increase in power and resolution afforded by whole genome approaches may challenge conclusions made on limited numbers of loci and maternally inherited haploid markers. Here, we provide a matched comparison of whole genome sequencing versus microsatellite and control region (CR) genotyping for Eurasian otters (Lutra lutra). Previous work identified four genetically differentiated “stronghold” populations of otter in Britain, derived from regional populations that survived the population crash of the 1950s–1980s. Using whole genome resequencing data from 45 samples from across the British stronghold populations, we confirmed some aspects of population structure derived from previous marker-driven studies. Importantly, we showed that genomic signals of the population crash bottlenecks matched evidence from otter population surveys. Unexpectedly, two strongly divergent mitochondrial lineages were identified that were undetectable using CR fragments, and otters in the east of England were genetically distinct and surprisingly variable. We hypothesize that this previously unsuspected variability may derive from past releases of Eurasian otters from other, non-British source populations in England around the time of the population bottleneck. Our work highlights that even reasonably well-studied species may harbor genetic surprises, if studied using modern high-throughput sequencing methods.


Genetic and Genomic Methods
As molecular ecology expands to include whole genome sequencing (WGS), the congruence between genetic and genomic methods has been called into question (McCartney-Melstad et al. 2018;Zimmerman et al. 2020).Rather than discrete categories, genetic and genomic methods form a spectrum across marker type, number of loci, and technologies used to generate data.Generally, methods using fewer than a thousand loci are considered genetic and therefore include microsatellites and mitochondrial fragments that are usually generated using Sanger sequencing-based methods (Hohenlohe et al. 2021).Methods using thousands or more loci are considered genomic, therefore including WGS and reduced representation methods such as restriction site-associated DNA sequencing (RADseq), which are produced using high-throughput sequencing technologies (Supple and Shapiro 2018;Hohenlohe et al. 2021).There is a trade-off in data acquisition between the number of individuals that can be assayed and the number of loci at which variation in each individual can be measured.WGS is becoming increasingly affordable but is still effectively limited to relatively small sample sizes.To maximize the potential benefit from the array of methods on the genetics-genomics continuum, an understanding of when these results are likely to be congruent, or differ, will enable the most costeffective marker selection for the specific research and conservation project (Gallego-García et al. 2021).
Ascertainment bias, the nonrandom analysis of loci resulting in parameter estimate biases (Nielsen 2000), has the potential to affect both genetic and genomic methods.For instance, microsatellite selection based on a single population will be biased toward detection of variation present in that population and away from detection of variation in a distinct, highly differentiated population (Malomane et al. 2018).This results in systematic Mol.Biol. Evol. 40(11):msad207 https://doi.org/10.1093/molbev/msad207Advance Access publication September 15, 2023 1 deviations, such as the underrepresentation of rare alleles, which can underestimate signals of a population expansion (Nielsen 2000).Genomic methods such as RADseq and WGS are less susceptible to ascertainment bias, due to the largely random choice of loci.Furthermore, if few loci show strong bias, these will have a larger effect in smaller, genetic data sets than in larger, genomic data sets.
Due to the varying evolutionary rates of different markers, more loci or greater length of sequence are required to obtain comparable levels of resolution and hence statistical power.For example, microsatellites are chosen to be multiallelic with an unusually high and variable mutation rate (from 10 −3 to 10 −4 mutations per locus per generation), whereas single nucleotide polymorphisms (SNPs) are primarily biallelic and have a lower mutation rate of 10 −8 -10 −9 per nucleotide per generation (Ellegren 2000; Väli et al. 2008).Microsatellites provide an effective method of identifying significant differences in variability even with low numbers of loci, and SNPs are often assessed in the thousands to millions in genomic data sets.Microsatellite loci are assumed to be largely neutral and thus useful for assessing genetic diversity and structure, whereas approaches such as WGS provide the opportunity to investigate the evolutionary mechanisms behind the measures of genetic variation, such as inbreeding depression, genetic basis of adaptation, and functional variation (Supple and Shapiro 2018;Hohenlohe et al. 2021).
It is important to evaluate characteristics of the study population when considering variation in results between methods.The variation found within species with large effective population sizes (e.g., N e > 1,000) is unlikely to be accurately quantified through analysis of a few loci.This issue is less critical in analyses of divergent, small, and inbred populations (Gallego-García et al. 2021).Furthermore, the number of samples taken to represent a population will constrain the power of analyses.For example, the potential to detect fine-scale population structuring is limited with few samples and the generally small sample sizes typical of genomic approaches may not yield the same insights as those made with much larger genetic sample size (Gallego-García et al. 2021).A clear benefit of genomic approaches is the ability to standardize results across research groups and therefore directly compare results between broad populations and species.For example, using a standard reference genome and short-read WGS data analyzed with the same bioinformatic pipeline facilitates direct comparison of data, which is rarely possible using microsatellite markers.
Case Study: Eurasian Otter (Lutra lutra) Eurasian otter (L.lutra) populations in Britain suffered a substantial population crash between 1950 and 1980 due to environmental chemical contaminants, such as persistent organic pollutants (Chanin and Jefferies 1978).This pattern was broadly observed across Europe, resulting in the species being classified as "near threatened" on the International Union for Conservation of Nature (IUCN) Red List (Roos et al. 2015).Although there are only limited direct estimates of otter population sizes throughout this time period, otter hunting records provide a proxy for population size.These indicate drastic population declines across southwest England and less drastic declines in northern England and Scotland from 1957 (Chanin and Jefferies 1978), highlighting a potential variation in the impact of chemical contaminants across the landscape.Restrictions and bans on the use of these chemicals in the 1970s lead to population recovery (Walker and Newton 1998;Mason and Macdonald 2004), as tracked by national surveys beginning in 1977 across Britain.The percentage of sites visited showing signs of otters increased from 6% to 59% in England, from 20% to 89% in Wales, and from 61% to 92% in Scotland, when comparing surveys from 1977 to 2009 for England and Wales and 1977 to 2003 for Scotland.The species has broadly been acclaimed as a conservation success story (Macdonald 1983;Strachan et al. 1990;Crawford and Scholey 2010;Kean and Chadwick 2021).However, the most recent National Survey of Wales identified signals of a recent, substantial population decline, highlighting the importance of continual monitoring and the value of otters as an indicator species for the ecological health of aquatic systems of Britain (Kean and Chadwick 2021).
Since the 1950-1980 population crash, extensive genetic assessments have been conducted on L. lutra in Britain, using both microsatellites (Dallas and Piertney 1998;Dallas et al. 1999Dallas et al. , 2002;;Hobbs et al. 2006Hobbs et al. , 2011;;Stanton et al. 2014;Thomas et al. 2022) and mitochondrial fragments (Stanton et al. 2009).These data have been incorporated into broader European assessments using both marker types (Cassens et al. 2000;Randi et al. 2003;Ferrando et al. 2004;Mucci et al. 2010).These studies illustrate the consequences of the population crash by revealing distinct genetic population strongholds.They were also critical in developing and sharing optimized microsatellite methods (Hobbs et al. 2011).A culmination of this work was a microsatellite data set spanning 21 years of the population recovery in Britain, which identified a time lag in the genetic connectivity between populations (Thomas et al. 2022).
Over 4,000 otters from across Britain, primarily victims of road traffic accidents, have been collected and preserved by the Cardiff University Otter Project since 1994, most of which are suitable for genomic analysis.A male L. lutra from Somerset, Southwest England, was sequenced and assembled to chromosomal completeness by the 25 genomes for 25 years project at the Wellcome Sanger Institute (Mead et al. 2020), and this high-quality reference genome is a strong backbone on which data can be mapped and contextualized.The Eurasian otter in Britain therefore presents an ideal opportunity to directly compare inferences derived from parallel data sets generated from genetic and genomic methods.As the first population genomic study of the species, we addressed two objectives.First, we used both WGS and an established microsatellite panel to directly compare population metrics derived from genomic and genetic methods for the du Plessis et al. • https://doi.org/10.1093/molbev/msad207MBE same sample set of Eurasian otters in Britain (fig.1).Secondly, we used the genomic data to assess inbreeding and historic effective population size estimates for the first time.We predicted that: 1) Detected genetic variation and population structuring will differ between genetic and genomic methods.2) The recent, anthropogenic population bottleneck (which was not previously identified using genetic methods) will be identifiable from genomic signatures, such as runs of homozygosity (ROH), linkage disequilibrium (LD), and changes in effective population size over time.

Results
Population We explored likely population substructuring in whole genome SNP and microsatellite data.The cross-validation (CV) error of ADMIXTURE analyses on the SNP data indicated the most likely value of K to be 1; however, the delta K of STRUCTURE analyses of the microsatellite data indicated the most likely value of K to be 3 (supplementary material S2.6, Supplementary Material online).Due to an increase in discordance among replicate STRUCTURE analyses and varying individual cluster assignments of the ADMIXTURE analyses at higher K values, results up to K = 6 are presented for both methods in fig.2d and e, whereas the results from the full range of K values are shown in supplementary material S2.6, Supplementary Material online.The STRUCTURE plot of microsatellite data indicated more admixture among populations than identified by ADMIXTURE on whole genome SNP data.Samples clustered broadly by geographic population in the STRUCTURE plot at K = 4, whereas in the ADMIXTURE analysis at K = 4, there was substructuring in the East, and Southwest England clustered with Wales.Southwest England and Wales were separated at higher K values.However, substructuring in the East was inconsistent from K = 4 to K = 6, with varying population assignments of the same individuals at increasing K values, consistent with the absence of pronounced substructuring within the East in the PCA (fig.2a).No pairwise relatedness of third degree or closer was identified among the samples (supplementary material S2.7, Supplementary Material online).
To probe the complex population assignments of individuals within the East population at increasing levels of K, fineSTRUCTURE analyses were run on the SNP data.The fineSTRUCTURE coancestry matrix (supplementary material S2.6, Supplementary Material online) showed clear population structuring between the East and the remaining populations, with less distinct population structuring between the Wales and North and Southwest England populations.The matrix also suggested lower levels of shared coancestry between the East and other populations, compared to within and between the non-Eastern populations.In contrast to the ADMIXTURE results, fineSTRUCTURE clustered mLutLut92 (East) with samples from the North.

Genomic Diversity within Britain
Based on whole genome SNP data, heterozygosity was higher in the East population than in other populations (mean difference of 0.00013-0.00017higher in the East; F 3,41 = 33.43,P < 0.0001; fig.3a).The realized inbreeding coefficient (F ROH ) based on the proportion of genomes in ROH ranged from 0.43 to 0.77 and was significantly lower in the East, with no significant differences among the remaining populations (F 3,41 = 16.21,P = 0.0004; fig.3a).The majority of ROH fell into the shortest, and therefore oldest class, indicating that extensive inbreeding occurred at 1,024 or more generations ago (an estimated 4,096 years).Signals of longer ROH (from 4 and 8 generations, or 16 and 32 years ago) were visible in a few individuals in the East and scarcer in the remaining populations.
Nucleotide diversity (π) was significantly higher across genomic windows in the East and significantly different among all population pairs except Southwest England and Wales (F 3,504 = 76.61,P < 0.0001; fig.3b).Private SNPs showed the same trend, with more SNPs identified as private to the East (2,076,325) than were private to any other population or common across all populations (1,498,476; fig. 3b; full private SNPs results in supplementary material S2.8, Supplementary Material online).
LD decay over distance (up to 100 kb) showed varying trends among populations (supplementary material S2.8, Supplementary Material online).Southwest England consistently showed the highest overall levels of LD and the slowest decline with increasing chromosomal distance; the North showed consistently lowest overall levels of LD, with the quickest decline over increasing chromosomal distance.Wales and East showed similar patterns of LD decay, which varied among chromosomes (supplementary material S2.8, Supplementary Material online).

Genetic Diversity within Britain
Genetic diversity statistics based on microsatellites were within the range of values reported in previous Genomics Reveals Complex Population History and Unexpected Diversity • https://doi.org/10.1093/molbev/msad207MBE microsatellite studies (Thomas et al. 2022) despite the lower sample size in the present study.Observed heterozygosity based on microsatellites was highest in the North and East (0.64 and 0.63), followed by Wales (0.53), and the lowest observed heterozygosity was identified in Southwest England (0.48).Since these statistics are not comparable across methods, we report the population rankings and full results in supplementary material S2.5, Supplementary Material online.When populations are ranked from highest to lowest heterozygosity, rankings varied between data sets: the whole genome SNP results rank East > North > Wales > Southwest England, whereas microsatellites rank North > East > Wales > Southwest England.Private alleles were identified in 12 of the 15 microsatellite loci, with 10 alleles across all loci private to the East, 11 alleles private to the North, and only 1 allele private to Southwest England.No alleles were found only in Wales (see supplementary material S2.5, Supplementary Material online).Rankings for F IS /F ROH showed a more notable difference, with whole genome SNP results ranking Southwest England > Wales > North > East and microsatellites ranking East > North > Southwest England > Wales.
Historic N e Historic effective population size (N e ) was estimated using GONE for the recent past (up to 800 years ago) and pairwise sequentially Markovian coalescent (PSMC) for the ancient past (10,000 to 1,000,000 years ago) (fig.4a; supplementary material S2.9, Supplementary Material Both the East and Southwest England populations showed substantial bottlenecks and recoveries between the 1950s and 1970s, but the N e of the Southwest was consistently higher than that of the East, which declined to 3.7 from Genomics Reveals Complex Population History and Unexpected Diversity • https://doi.org/10.1093/molbev/msad207MBE 1972 to 1984 (averaged across bootstrap N e estimates for 9-12 generations ago).The decline in N e started earlier in Wales, during the 1800s, and showed an increase in the past 50 years.The North showed a gradual, continual decline through the past 800 years.For the most recent 50 generations, we compared N e estimates using GONE (fig.4b) to the survey data recording the number of sites showing positive signs of otters by region (fig.4c).Visually, the trends matched reasonably well, with the differences in trends between regions and the extremely low proportion presence in the East coinciding with the low N e estimates.At deeper timescales, PSMC analyses of 10,000 to 1,000,000 years ago showed a decline in N e across all populations, albeit with some local variation.

Mitochondrial Genomes versus Control Region
Of the 45 British samples, whole mitochondrial genomes for 44 were assembled, of which 41 were assembled using NOVOPlasty and 3 were assembled using MITObim (table 1).From these 44 sequences, 18 unique haplotypes

MBE
were identified across 153 segregating sites (summary statistics and haplotypes given in supplementary materials S1 and S2.11, Supplementary Material online, respectively).The TCS network identified 2 distinct lineages, equivalent to lineages 1 and 3 in the lineage classification of Waku et al. (2016), separated by a branch representing 105 mutations between groups (fig.5a).Due to the presence of both divergent lineages in the East, the number of segregating sites and nucleotide diversity were higher in the East than other British populations; however, haplotype diversity was highest in the North (supplementary material S2.11, Supplementary Material online).
All five control region (CR) haplotypes which have previously been identified in Britain were found in our data set.Haplotypes Lut1, Lut3, Lut6 (Stanton et al. 2009), Lut4 (Cassens et al. 2000), and Lut7 (Pountney 2008) Genomics Reveals Complex Population History and Unexpected Diversity • https://doi.org/10.1093/molbev/msad207MBE were found in 29, 4, 9, 1, and 1 samples, respectively.The difference between Lut1 and Lut4 is a single base indel at nucleotide position 96, and therefore, these haplotypes were collapsed in the haplotype network (fig.5b).Geographic distribution of these haplotypes was broadly consistent with previous findings.For example, Lut1 was found across all populations, whereas Lut3 was only identified in the East.Interestingly, however, Lut6, previously only found in Western Britain, was identified in eight samples from Wales and one sample from the East in our data set.Only the four samples belonging to lineage 1 corresponded to CR haplotype Lut3, separated from the remaining haplotypes by a single mutation.
We also assembled mitochondrial genomes from three previously generated short-read data for European L. lutra available in the short-read archive.These read sets were successfully assembled using NOVOPlasty and aligned with European and Asian L. lutra mitochondrial genome sequences from GenBank (n = 13) and the British samples generated in this work, yielding a total of 60 L. lutra sequences.We added a single hairy-nosed otter (Lutra sumatrana) as outgroup.From the 60 L. lutra mitochondrial genomes, 34 unique haplotypes were identified across 772 segregating sites.Phylogenetic analysis showed separation of samples into four main lineages, all with posterior probabilities of 1.00 (fig.5d; sample lineages in supplementary material S2.11, Supplementary Material online).Three lineages were named by Waku et al. (2016) as L. l. nippon (here referred to as lineage 5), lineage 1, and lineage 2, and one new lineage was identified in this study (lineage 3).Geographic locations of these lineages are provided in fig.5c, with the exception of Britain and Russia for which multiple lineages were identified.The British samples sequenced in this study identified two divergent lineages: lineage 1 and lineage 3. Sequences assigned to lineage 3 include 40 samples from Britain (across all populations) and the SRA-derived sequences from Denmark, Russia, and Norway.Lineage 1 included four samples from the East of Britain, alongside GenBank samples from China, Laos, and Japan, and the reference genome sequence, from a British otter (from Somerset, closest geographically to our Southwest England region).Although the rooted phylogenetic tree indicated that lineage 5/L.l. nippon split first, followed by lineage 3 and then lineages 1 and 2, variation in sequences within lineage 3 appeared to be a recent diversification relative to the older branching within lineages 1 and 2. These results were reflected in the higher number of segregating sites and nucleotide diversity observed within lineages 1 and 2, relative to lineage 3.

Genetic versus Genomic Methods
Our study provides a direct comparison of genetic and genomic methods in assessment of genetic diversity and population structure in a threatened wild carnivore, the Eurasian otter.We applied genomic methods to investigate the occurrence and effects of a recent, anthropogenically driven population decline and recovery of Eurasian otters in Britain.Our concurrent analysis of microsatellite and SNPs from the same sample set highlights the complexities inherent in interpreting results of these approaches.
Broadly, our microsatellite and SNP data sets are not concordant, with substantial differences in the order of population differentiation and population groupings.For example, at K = 3, the SNP data set grouped Southwest England and Wales, whereas the microsatellite data grouped East and North England.Across all samples, the proportion of admixture identified is higher in the microsatellite data set than the SNP data.However, for the

MBE
purpose of population assignment based on a suitable K value, there is also some concordance across data sets, with the exception of some samples between the East and North in the microsatellite data.To summarize, similarities in the ADMIXTURE/STRUCTURE comparison broke down at higher values of K and in more complex scenarios such as the substructuring within the East.
Similarly, SNPs showed more population structuring in PCA, and microsatellites showed significantly higher pairwise F ST estimates.These results are not unexpected, as both similarities (Zimmerman et al. 2020) and differences (Lah et al. 2016) such as these have been identified in past studies of comparable genetic and genomic data sets.
Both the number of loci assessed and ascertainment bias are likely to be contributing to the differences between the genetic and genomic approaches.Specifically, the higher resolution of data captured by almost 9 million Genomics Reveals Complex Population History and Unexpected Diversity • https://doi.org/10.1093/molbev/msad207MBE biallelic SNPs, compared with 15 multiallelic microsatellites, is likely to be critical, as in several prior studies (Lah et al. 2016;Natesh et al. 2017;Lavretsky et al. 2019;Gallego-García et al. 2021).Although the microsatellites identified the broad patterns of population structuring accurately, they do not hold the same power to detect more fine-grained population distinction (higher values of K) and more complex scenarios, when compared with the genomic results.Since the complexity of any given study system is typically unknown, it is important to recognize that conclusions based on genetic methods may change following the application of genomic data.In particular, the measures which showed the most significant differences between methods (e.g., pairwise F ST ) are likely to be inflated due to ascertainment bias of the microsatellite loci.Specifically, these microsatellite loci were identified and selected based on the variation they identified in Scotland and Wales (both sampled in this study), and therefore, they are likely to represent only a subset of the variation present across populations in Britain (Dallas and Piertney 1998;Dallas et al. 1999;Hobbs et al. 2006).Although ascertainment bias is less commonly identified as a source of variation between genetic and genomic comparisons, Fischer et al. (2017) showed significantly larger estimates of pairwise F ST using microsatellites than SNPs, which they attributed to the ability of high-coverage WGS in identifying rare variants, which reduce overall F ST .Similarly, Cairns et al. (2023) identified higher proportions of admixture among dogs and dingoes (Canis spp.) using microsatellites relative to SNPs, where the microsatellites had been selected based on allele frequencies identified in dogs.Taken together, these findings highlight potential impacts of microsatellite ascertainment bias, leading to inflated admixture proportions (Cairns et al. 2023).
Our utilized microsatellite panel does not contain any loci pairs with significant LD signals, and all loci are at least 1.34 Mbp apart in the genome.Consequently, linkage appears unlikely to have had a large effect on our microsatellite results.Linkage signals are higher for our analyzed SNPs, due to their larger numbers (9 million) and greater proximity in the genome.Therefore, although linkage occurs in both data sets, it is likely to be affecting the genomic results more significantly than the genetic results in this study.At this point, it becomes difficult to disentangle the effects of number of loci, ascertainment bias, and linkage, as using too few loci or using nonrandomly selected loci could both underrepresent rare alleles and therefore lead to inflated F ST estimates, among other effects, further supporting the crucial role of genomic data.
We are not suggesting microsatellite markers are redundant, and we emphasize the important role they continue to play as an affordable marker system for cluster assignment and other scenarios.Equally, long-term genetic data sets, such as Thomas et al. (2022), highlight the importance of the continuity of markers to assess temporal trends in genetic variation, for example, during population recovery and expansion.However, our analyses suggest that due to numbers of loci and ascertainment bias, the absolute values of metrics calculated based on microsatellites should be interpreted with caution, with a focus on relative values (such as a decrease in F ST over time) likely being much more informative.Importantly, the scale of the effect of ascertainment bias varies depending on the specific loci and samples, and therefore, comparisons between studies based on different loci, populations, samples, or species would be inappropriate.
Using whole mitochondrial genomes, we were able to identify two divergent lineages within Britain, which were not previously identified based on CR studies (Stanton et al. 2009), and were not possible to identify using CR fragments for our samples.Based on these results, we strongly recommend that analyses using short mitochondrial fragments should be upgraded to whole mitochondrial genome analyses where possible, to avoid misleading inferences.

Demographic Analyses and Consequences
Using whole genome sequence data, it was possible to detect and analyze historic trends in effective population size and access methods not available for genetic data sets.Most notably, the recent anthropogenic population bottleneck between the 1950s and 1970s was clearly evidenced (using GONE), with a sharp decrease and subsequent increase in N e occurring concurrently across southern England (based on samples distributed across southwest, central, and eastern England).In contrast, Wales showed a decline beginning around 1830s, coinciding with the "improvement of the sporting gun with percussion detonation," and reduced otter hunt results in North Wales (Jefferies 1989).N e estimates from Wales increased from the 1950s (as other British populations were beginning to decline), which was unexpected but potentially due to the impact and changes of otter hunting around this time, despite legal protection only being obtained in 1978 in Wales and England (Jefferies 1989).The North showed an overall, gradual decline from around 1200 to the present day, with no change in this trend through the 1950s-1970s bottleneck detected in southern England, confirming suspicions that this larger, more rural population was less affected by chemical pollution (Chanin and Jefferies 1978).These findings are broadly in keeping with otter hunting records which, for example, indicate a decline in success rate of otter hunts from 1957 across southwest England, a lesser decline in northern England and Scotland, and no evidence of a decline in north Wales (Chanin and Jefferies 1978).The accompanying PSMC analyses also show a more long-term decline from 1,000,000 to 10,000 years ago across all populations.
We hypothesized that genomic signatures of a population bottleneck would be evident.Breeding between relatives leads to homozygous regions of the genome separately inherited from a recent shared ancestor, known as ROH.Longer ROH indicate recent inbreeding, and shorter runs indicate older inbreeding, which have been broken up by recombination (Druet and Gautier 2017).We predicted that the British bottleneck occurred du Plessis et al. • https://doi.org/10.1093/molbev/msad207MBE between the 1950s and 1970s and would therefore be identifiable as ROH of a length corresponding to inbreeding occurring 9-17 generations before these samples were collected (based on a generation time of 4 years); however, ROH of this length were infrequent.ROH have been identified in populations suffering severe reductions in effective population size for prolonged periods of time, such as the Florida panther, Puma concolor (Saremi et al. 2019).We do not have census data spanning the bottleneck, but due to the banning of some chemical contaminants in the 1970s (Walker and Newton 1998), we predicted that the bottleneck was unlikely to have persisted for more than 20 years.However, the GONE analyses indicated that depending on the population, N e dropped for between 30 and 45 years before increasing (with the exception of Wales and the North).Despite the length and severity of this bottleneck, especially in the East (from the lowest N e estimate of around 100 to around 3), it does not appear to have led to a consistent burden of ROH across multiple individuals.Conversely, we observe many, short ROH across all samples.Theory would suggest this indicates very old inbreeding followed by generations of recombination, potentially from a bottleneck that is older than we are able to estimate using RzooROH or due to background relatedness (Ceballos et al. 2018).We found no evidence for bottlenecks from 10,000 to 1,000,000 years ago (PSMC analyses) nor from 200 to 800 years ago (GONE analyses); however, we have not assessed evidence for the period 800 to 10,000 years ago due to method limitations.It is unexpected, but reassuring, that we do not see long ROH in these individuals despite the severe recent bottleneck in the East.However, the evidence for extensive historic inbreeding and its legacy, the generally high realized inbreeding coefficient (F ROH ) observed in the modern populations, provides cause for concern with respect to the species' long-term viability, particularly given the current context of small populations and anthropogenic threats for a near threatened species (Roos et al. 2015;Reid et al. 2019).

Distinct Signatures Detected in the Eastern Population
Unexpected diversity identified in the East using genomic methods was not previously identified in studies using genetic methods (Stanton et al. 2014).Estimates of variation (e.g., from PCA) were dominated by variation among samples within the East, and these results were matched by higher heterozygosity, nucleotide diversity, and private SNP counts, alongside lower inbreeding coefficients relative to other British populations.The GONE analyses indicate that the population from the East went through the most severe bottleneck of all British populations, with N e as low as 3.7, making the high variation in this region particularly surprising.Although the ROH analyses identified some recent inbreeding, it was not consistent across individuals from the East, indicating a more complex demographic history than hypothesized from a simple anthropogenic bottleneck.Other metrics are also consistent with a relatively complex demographic history in the East population, including ADMIXTURE results, which indicate inconsistent population clustering, and fineSTRUCTURE results, which indicate that the East have lower shared coancestry with remaining populations than they show when compared with one another but also show inconsistent coancestry between individuals from the East.For example, although some pairs within the East show very high shared coancestry, others show very low shared coancestry, whereas all individuals in Wales show a very similar shared coancestry between pairs.Lastly, the divergent mitochondrial lineage identified using genomic methods was found only in samples from the East and the reference genome (from Somerset, equidistant to southwest and eastern samples used in this study).In line with previous studies (Thomas et al. 2022), SNP-based measures of pairwise F ST were not larger for pairs including the East, suggesting that although the East contains lots of unique variation, it is not more distantly related to remaining populations.Since pairwise F ST is a relative measure of genetic variation within a population relative to total variation, it is possible that these results reflect the high variation observed within the East, when compared with total variation.
The results observed in the East do not align with the simple scenario of a population bottleneck and recovery, and instead, we propose alternative hypotheses to explain our findings.In the first National Otter Survey of England, in 1977England, in -1979, of 623 , of 623 sites surveyed in East Anglia, only 20 were positive, illustrating how close the population was to local extinction (Lenton et al. 1980).Captive-bred otters were released, 13 in the south and 81 in the east of England, as a part of a broader reintroduction program of over 180 otters across Britain (Green 1997).Following these releases, the second National Survey in 1984-1986 found 8 positive out of 725 surveyed sites in East Anglia, 5 of which were from released otters (Strachan et al. 1990).One explanation for genetic distinctiveness of otters in the East might be that their small effective population size led to strong genetic drift.However, the high genetic diversity (both overall variability and the high private SNP count) observed in the genomic data in the East renders this explanation unlikely.Another explanation is that the East gene pool contains ancestry from Eurasian otters of non-British origin, explaining both the unusually high proportion of unique variation and the divergent mitochondrial lineage.A study by the European Association of Zoos and Aquaria (EAZA) using microsatellites identified two main genetic "lines" of Eurasian otters in captivity, known as A-and B-lines (E.Rey pers.comm.),where anecdotal evidence suggests the B-line otters were bred with the Asian subspecies L. l. barang found in Sumatra, Thailand, and Vietnam (J.Palmer pers.comm.(Hung and Law 2016)).An organization involved in the reintroduction program was separately breeding Eurasian otters in captivity using two founders from Thailand, thought at the time to be L. l. barang (J.Palmer pers.comm.).Although studbook records show no evidence of crossing between otters of Thai origin and British L. lutra (J.Palmer pers.comm.),Genomics Reveals Complex Population History and Unexpected Diversity • https://doi.org/10.1093/molbev/msad207MBE escapes or unknown mixing between pens cannot be excluded, however unlikely (P.Chanin pers.comm.).Although there are no whole mitochondrial genome sequences available from Eurasian otters from Thailand, we have included a sequence from the neighboring Laos in this study, which groups with samples from China and the east of England.Therefore, when combined with our genomic results, these findings indicate the possibility of individuals of Asian origin (likely Thailand) being either accidentally bred, released, or escaping into the south and east of England.Otter releases in this region occurred between 1983 and 1996, coinciding with the increase in our estimate of effective population size for the East.We therefore suggest that the introduction of a few Asian-origin or admixed otters to a very small existing population had a large impact on the genetic identity of the population, leading to the high proportion of unique genomic variation and divergent mitochondrial lineages.The captive B-line otters are no longer being bred in captivity (J.Palmer pers.comm.);however, there are indications that some of their descendants have been released in Europe (Hájková et al. 2007).Our results show that the consequences of such introductions have left genomic signatures across the east of England and possibly beyond.
To provide further evidence describing the genetic history of otters in Britain, we suggest two main paths of future work.Firstly, to further investigate the history of otters in the East of Britain, it would be beneficial to sequence historical samples from the region (taken before the population bottleneck), as well as B-line captive-bred otters, to compare these to the population we observe presently and assess the likelihood of a replacement event occurring through the reintroductions (Strachan et al. 1990).Secondly, we recommend sampling the range of the Eurasian otter broadly, including historic samples, to try to identify any similarities between a potential source population for British reintroductions, such as L. l. barang, L. l. chinensis (found in Thailand), and other Asian lineages, alongside any other populations where B-line otters may have been reintroduced.This would also enable a thorough genomic assessment of subspecies classifications, clarifying the existing confusion around Asian subspecies, such as L. l. barang and L. l. chinensis.We note that population declines and local extinctions across the species' range mean that it is possible that the source population could be missed in this sampling, highlighting the importance of historic samples (Yoxon and Yoxon 2019).

Conservation Implications
Our genetic to genomic comparison highlights the value of genomic methods in conservation to avoid misinterpretations of potentially biased, low-resolution markers.However, this does not make the Eurasian otter microsatellite panel or previous studies based on this redundant.Rather, this study identifies the limits of the interpretations of these results when compared with genomic data.A SNP array for capture-based SNP genotyping could be designed based on our findings, to provide cheaper genomic analyses of British otters.However, based on our findings of ascertainment bias, we would recommend a SNP array based on SNPs identified from a range-wide sampling of the species.
Here, we provide the first clear evidence of a population bottleneck in Eurasian otters in southern England, highlighting an important contrast to the history of populations in northern England, Scotland, and Wales.We provide clear evidence of previously unexpected and unusual signatures in the east of England and assign these to a mitochondrial lineage only found, as yet, in Asia.Both of these findings have only been possible through the application of genomic methods, enabled by the prior publication of an extremely high-quality reference genome (Mead et al. 2020), illustrating the importance of WGS for conservation of this and other species.

FIG. 1 .
FIG. 1. Geographic origins of Eurasian otter samples analyzed in this study and the reference genome.

FIG. 2 .
FIG. 2. Comparison of estimates of population structuring of Eurasian otters in Britain, based on WGS versus microsatellite data.PCAs of (a) whole genome SNPs and (b) microsatellites for the whole data set and (in the inset) for the 33 non-Eastern samples and (c) pairwise F ST values for whole genome SNPs and microsatellites.(d ) ADMIXTURE results for whole genome SNPs and (e) STRUCTURE results for microsatellites.

FIG. 3 .
FIG. 3. Genomic diversity in British Eurasian otters based on SNP markers.Individual-level (a) genome-wide heterozygosity and proportion of the genome in runs of homozygosity (F ROH ) as calculated by RzooROH, color-coded by number of generations ago inbreeding is estimated to have occurred (from older runs/shorter ROH to recent/longer ROH).Population-level measures (b): nucleotide diversity (π) calculated in 20 Mb sliding windows and counts of SNPs private to each population (bars), relative to the number of SNPs common to all populations (black line at 1,498,476 SNPs).

FIG. 4 .
FIG. 4. Estimated effective population size through time by population of Eurasian otter in Britain based on SNP data and comparable survey data.(a) GONE analyses were conducted on a population-level scale, with ten bootstrap repeats for each population, assuming a generation time of 4 years, conducted on samples sizes of n = 8 for Southwest England, n = 12 for Wales and East, and n = 13 for North, and the expected timing of the bottleneck highlighted.(b) GONE analyses for the same period of time for which surveys were conducted (0-50 years before samples were collected in 2020).(c) Survey data as percentage of sites surveyed which showed positive signs of otters, grouped by regions reflecting stronghold populations, with the exception of the North, which was split into regions from England and Scotland based on divergent census size trends.For surveys which spanned more than 1 year, the results are plotted at the earliest year.All plots are presented on a log scale of N e .

FIG. 5 .
FIG. 5. TCS networks, phylogeny, and locations of lineages of mitochondrial genome variation in Eurasian otters.(a) TCS network of 44 British samples and the reference mitochondrial genome (LR822067.1),based on 16,365 bp mitochondrial sequence (with the repeat region removed).(b) TCS network of the same 44 British samples and known CR haplotypes from GenBank (in black).Lut4 and mLutLut29 are collapsed into Lut1, as they only differ by a single base indel.(c) Countries shaded to indicate where each whole mitochondrial genome lineage was identified, and hatching indicates known, global Eurasian otter range.Japan included both lineages 1 and 5 and therefore is shaded orange for the unique lineage 5.No colored shading is shown for Russia and Britain in this diagram, due to multiple lineages identified (see supplementary material S2.11, Supplementary Material online).(d ) Phylogenetic tree of 16,392 bp whole mitochondrial sequence (with the repeat region removed) generated in this study (n = 44), assembled from SRA (n = 3), and downloaded from GenBank (n = 13), rooted with a hairy-nosed otter (L.sumatrana), totaling 61 sequences.Posterior probabilities of branches between lineages are all 1.Tree shading indicates lineage 1 (red), 2 (green), 3 (blue), and 5/L.l. nippon (orange).Full phylogeny and posterior probabilities are given in supplementary material S2.11, Supplementary Material online.

Table 1 .
Eurasian Otter Sample Accessions, Locations, and Publication Source of Publicly Available SRA and Whole Mitochondrial Genome Sequences Incorporated into Analyses.