Mitochondrial genome-based synthesis and timeline of Eurasian otter (Lutra lutra) phylogeography

ABSTRACT

Eurasian otters (Lutra lutra) have a broad distribution across Eurasia, but biogeographic data outside of western Europe is limited to disconnected pockets. Based on current subspecies designations, Asia appears to harbour a large proportion of the species’ diversity, with 10 of 12 Eurasian otter subspecies found in Asia. Here we provide a range-wide synthesis of mitochondrial data, inferring a timeline and pattern of phylogeographic signals. Whole mitochondrial genomes of 27 Eurasian otters across 4 subspecies are presented from newly generated data (n = 6; 4 from Korea, 1 from Hong Kong and 2 from UK), assembled from the Sequence Read Archive (n = 4), and sourced from GenBank (n = 17). We then combined whole mitochondrial genome results with cytochrome b data to increase the sample size and contextualise our results with prior studies. We identified five distinct lineages that were discordant with current subspecies classification. Phylogenetic dating revealed that the earliest diverging lineage was the Japanese lineage, with remaining lineages diverging ≥0.08 million years ago. Mitochondrial diversity calculated by sample locations seemed mainly driven by the presence of multiple lineages. When grouping samples by lineage, genetic diversity was highest in Lineage 1 (primarily found in China and Laos), followed by Lineage 2 (primarily found in Korea), and lowest diversity identified in Lineage 3 (primarily found in Europe). Our findings highlight previously undetected lineage diversity within Eurasian otters, but also the need for further taxonomic and genomic evaluation of the species in Asia. The identified unique, distinct lineages of Eurasian otters also warrant urgent conservation attention.

KEYWORDS:

• Biogeography
• cytochrome b
• mitochondrial genome
• phylogeography
• subspecies

Introduction

Eurasian otters (Lutra lutra; Linnaeus, 1758) have the largest geographic distribution of any otter (subfamily Lutrinae) species, spanning from the Iberian Peninsula in the West, to Japan in the East, from Norway and Russia in the north, to North Africa and Southeast Asia in the south (Roos et al. 2015). Population declines and bottlenecks have prompted protective legislation and extensive monitoring and research in Western Europe, but only limited protection or research has been extended to subspecies and populations in Asia (Basnet et al. 2020), despite reports of severe declines in China during recent decades (Li and Chan 2018; Chen et al. 2023). Threats to this globally ‘near threatened’ species vary across their range and are broad, encompassing chemical pollution, habitat loss and fragmentation, and illegal trade (Roos et al. 2015). As in other widely distributed Eurasian carnivores such as brown bears (Ursus arctos) and red foxes (Vulpes vulpes) (Davison et al. 2011; Statham et al. 2014), it is likely that these understudied regions of Asia hold a large proportion of the diversity of the species, both genetically and phenotypically. Indeed, 10 of the 12 known Eurasian otter subspecies are found in Asia (Hung and Law 2016).

Despite the debate around delineating populations into evolutionarily significant units, subspecies and species, these classifications and terms have value in the practical conservation and management of endangered species (Fraser and Bernatchez 2001). Multiple subspecies of Eurasian otters have been described over the last 250 years based on a combination of geographic distribution and morphological evidence, usually derived from a limited number of examined specimens. Between seven and 28 subspecies have been cited, with the 2015 IUCN Red List report indicating seven subspecies (Roos et al. 2015), although another commonly referenced paper indicates 12 extant subspecies (Hung and Law 2016). The only genomic assessment of Eurasian otter subspecies was conducted by Waku et al. (2016), identifying the genetically distinct Japanese otter using whole mitochondrial genome sequences. Currently, the data required to establish the Japanese otter as a separate species (Lutra nippon) or a Eurasian otter subspecies (Lutra lutra nippon) is lacking, for example phenotyping and morphology, and therefore we refer to the Japanese otter as a subspecies, L. l. nippon as in Waku et al. (2016). Publication of the evidence for the distinct Japanese otter in 2016, four years after Eurasian otters were declared extinct in Japan in 2012 (Waku et al. 2016), highlights the urgency of further genetic assessments of the species. Many Eurasian otter subspecies are currently restricted to small geographic regions, for which only very limited information on the species’ conservation status is available. For example, L. l. aurobrunnea is endemic to ‘the lower and central hilly region of Nepal’ (Hung and Law 2016), where there have been ‘no recent sightings documented’ and ‘literature reports [are] largely over 30 years old’ (Yoxon and Yoxon 2019). Although a combination of risk factors is leading to global otter population declines, the structuring of small, remote, isolated populations enhances the risk of putative subspecies, such as L. l. aurobrunnea, going extinct before any research on the presence or taxonomy of the subspecies is conducted.

Based on a lack of sightings since 1979, L. l. nippon was declared extinct in Japan in 2012. However, Eurasian otters were recorded on camera traps in Japan in 2017 (Nakanishi and Izawa 2019). Nakanishi and Izawa (2019) hypothesise that otters may have naturally or artificially been reintroduced to Japan across the 60 km strait from the coastal islands of the Korean Peninsula. If otters did manage to cross from Korea to Japan, gene flow of L. l. nippon from Japan to Korea – before the Japanese population had gone extinct – is also plausible. However, prior studies of mitochondrial variation in Korea have shown very low diversity (Koh et al. 2004), indicating that this region is unlikely to harbour any L. l. nippon. As one of Japan’s nearest neighbouring countries, analysis of multiple samples from Korea provides the opportunity to investigate within-region variation, while also evaluating the possibility that the Japanese lineage (L. l. nippon), currently believed to be extinct, might in fact persist outside of Japan.

Genetic tools provide an insight into the evolutionary history, processes, and divergences among populations, offering a framework for assigning subspecies classifications based on phylogenetic history. Commonly used markers include control region and cytochrome b regions of the mitochondrial genome, often ranging from 100–400 base pairs (bp) in length (Cassens et al. 2000; Ferrando et al. 2004), with some examples of fragments 900–1600 bp in length (Mucci et al. 1999; Koh et al. 2004; Koepfli et al. 2008; Coudrat et al. 2022). In Eurasian otters, control region fragments have been utilised most extensively in previous studies, however this fragment has a very limited resolution (du Plessis et al. 2023). In contrast, cytochrome b has been less frequently used, but as a 1140 bp fragment, it potentially has a higher resolution for fine-scale, recent evolutionary processes. This marker has been used in phylogenetics (Koepfli and Wayne 1998) and species identification (Fernandes et al. 2008), but more recently has been effective in assessing regional taxonomy and phylogeography, identifying haplotypes from Italy, Iraq (Moretti et al. 2017), Korea, Japan (Park et al. 2019), and Laos (Coudrat et al. 2022). Even fewer studies have used whole mitochondrial genome sequences, at least 16,000 bp in length (Jang et al. 2009; Waku et al. 2016), which provide a significant improvement in genetic resolution (i.e. number of variable sites, and hence statistical support for phylogenetic groupings; (Keis et al. 2013)) over both control region and cytochrome b fragments. Increasingly, studies are using concatenations of multiple mitochondrial markers, resulting in 1500–4000 bp in length, analysed in parallel to autosomal microsatellites or a subset of whole mitochondrial genomes (Coudrat et al. 2022; Jang-Liaw et al. 2023). Both Coudrat et al. (2022) and Jang-Liaw et al. (2023) highlighted the value of local studies, using established markers in combination with a subset of whole mitochondrial genomes. Specifically, they have been applied to species identification and occurrence, and baseline population genetics in regions where conducting Eurasian otter population genetics is challenging and existing knowledge is currently limited (Coudrat et al. 2022; Jang-Liaw et al. 2023).

Here we present a phylogeographic assessment of the Eurasian otter, based on whole mitochondrial sequences across the species range, with a focus on populations in Asia. To contextualise our findings with prior studies, we also include analyses of the cytochrome b gene, which in previous studies has been surveyed across a broader geographic range. Combining evidence from high-resolution whole mitochondrial genome sequences with the larger sample size of cytochrome b sequences available, we address the following three hypotheses: (1) that the distinct Japanese lineage (L. l. nippon) is not present in Korea or elsewhere in eastern Eurasia; (2) that genetic lineages do not reflect subspecies classifications based on geography; and (3) that there is higher mitochondrial variation in Asia than in western Europe.

Methods

Sample collection and sequencing

Samples (Table 1) from UK (GBR01, GBR02), Hong Kong (CHN01) and South Korea (KOR01, KOR02, KOR03) were collected and sequenced using short, paired end reads, for this study (see Supplementary Material).

Table 1. Eurasian otter sample accessions, locations, and publications for whole mitochondrial genome data included in this study. Metadata for all newly generated, publicly available SRA and whole mitochondrial genome sequences analysed in the study. Whole mitochondrial genome lineages are shown, as inferred in this and prior studies.

Sample accession/ID	Location	Publication	Subspecies	Lineage
Newly sequenced (this study):
CHN01	Hong Kong, China	This study	L. L. chinensis	1^+
KOR01	South Korea	This study	L. l. lutra	2^+
KOR02	South Korea	This study	L. l. lutra	2^+
KOR03	South Korea	This study	L. l. lutra	2^+
GBR01	Northeast England, UK	This study	L. l. lutra	3^+
GBR02	Central England, UK	This study	L. l. lutra	3^+
SRA reads downloaded and assembled:
SRR19383068	Narvik, Norway	de Ferran et al. (2022)	L. l. lutra	3
SRR19383067	Tyumen Oblast, Russia	de Ferran et al. (2022)	L. l. lutra	3
SRR11679564	Denmark	Margaryan et al. (2021)	L. l. lutra	3
SRR5192239	Sumatra	Unpublished^1	L. l. barang	4^+
Whole mitochondrial genome sequences from Genbank:
LC049377	China	Waku et al. (2016)	–	1
LC049378	China	Waku et al. (2016)	–	1
LC049952	Sichuan, China	Waku et al. (2016)	L. l. chinensis	1
LC049953	Unknown	Waku et al. (2016)	–	1
LC049954	Sakhalin, Russia	Waku et al. (2016)	L. l. lutra	2
LC049955	Kanagawa, Japan	Waku et al. (2016)	L. l. nippon^2	1
LC050126	Kochi, Japan	Waku et al. (2016)	L. l. nippon^2	5^+
LC094961	Laos	Waku et al. (unpublished)	L. l. barang	1
LR822067	SWEng, UK	Mead et al. (2020)	L. l. lutra	1
MW316682	Kinmen Island, Taiwan, China	Jang-Liaw et al. (2023)	L. l. chinensis	1
EF672696	Korea	Ki et al. (2010)	L. l. lutra	2
FJ236015	Korea	Jang et al. (2009)	L. l. lutra	2
MW573979	Daejeon, South Korea	Kim and Jo (2021)	L. l. lutra	2
MW344881	Fuzhou, China	Jang-Liaw et al. (2023)	L. l. chinensis	1^+
OP554563	Nakai-Nam Theun National Park, Laos	Coudrat et al. (2022)	L. l. barang	1^+
OP554564	Nakai-Nam Theun National Park, Laos	Coudrat et al. (2022)	L. l. barang	1^+
OP554565	Nakai-Nam Theun National Park, Laos	Coudrat et al. (2022)	L. l. barang	1^+
KY117556	–	Mohd Salleh et al. (2017)	–	–

¹ Raw reads published accompanying Mohd Salleh et al. (2017), but mitochondrial genome assembly was unsuccessful in the original study.

² Classified as a separate species, Lutra nippon, or subspecies, Lutra lutra nippon (Waku et al. 2016), and by Gray (1867) referred to as the subspecies, Lutra lutra whiteleyi.

⁺ Lineages inferred in this study.

Previously published whole mitochondrial genome data

Paired-end Sequence Read Archive (SRA) data was available for three samples, from Russia, Norway and Denmark, alongside single-end read data for one sample collected in 1923 from Sumatra (Table 1). These data were downloaded and assembled using the same pipeline as the newly generated data, except for the Sumatran sample (SRR5192239). A further 16 previously published whole mitochondrial genome sequences were available on GenBank alongside the one from the Eurasian otter reference genome (LR822067), and the Hairy-nosed otter (Lutra sumatrana) was included as an outgroup, leading to a total of 28 whole mitochondrial genome sequences included in these analyses (Table 1). Known sample locations and mapping for regional samples are given in Supplementary Material.

Putative subspecies classifications were assigned based on the known locations of samples and the geographic distribution of subspecies according to Hung and Law (2016), where possible. Currently, there is no other information which could be used to discern subspecies. From the 12 currently recognised subspecies, we have sampled four; L. l. lutra, L. l. chinensis, L. l. barang and L. l. nippon.

Whole mitochondrial genome assembly and analyses

All newly generated and SRA data were assembled using NOVOPlasty v4.3.1 (Dierckxsens et al. 2017), except for sample SRR5192239 which was assembled using MITObim v1.9.1 (Hahn et al. 2013) (see Supplementary Material).

All whole mitochondrial genome sequences were aligned using MAFFT v7.481 (Katoh et al. 2005; Katoh et al. 2002), using the ‘ – auto’ option to automatically select an appropriate alignment strategy, which was FFT-NS-2 (fast, progressive method). Sequences were aligned with and without the outgroup, before removing the repeat region in Geneious Prime v2022.0.2 (https://www.geneious.com). Both with and without the outgroup aligned, the repetitive region was removed from bases 16,034–16,255 (relative to the reference mitochondrial sequence, LR822067), resulting in a final alignment length of 16,471 bp with the outgroup, and 16,432 bp without the outgroup. PopArt v1.7 (Leigh and Bryant 2015) was used to produce a TCS network (Clement et al. 2002) of this alignment, with the outgroup removed, and indels not considered (therefore some haplotypes are collapsed into a single node). The R packages pegas (Paradis 2010) and ape (Paradis and Schliep 2019) were used to calculate summary statistics including segregating sites, haplotype diversity and nucleotide diversity (π). The standard error of the nucleotide diversity estimate was calculated as the square root of the variance estimated by pegas, divided by the square root of sample size. Summary statistics were calculated by groups based on both broad geographic region and whole mitochondrial lineage membership (see SM1). BEAST2 v2.6.6 (Bouckaert et al. 2019) was used to produce a phylogeny, and date divergences. ModelFinder (Kalyaanamoorthy et al. 2017) within IQ-TREE (Nguyen et al. 2015) was used to identify the best-fitting model, based on the Bayesian Information Criterion (BIC), which could be implemented in BEAST2. This was the HKY model (Hasegawa et al. 1985), defined as unequal transition/transversion rates, using unequal, empirical base frequencies, and estimating the proportion of invariant sites. Variation across sites was modelled using a gamma distribution, with a kappa (ratio of transition and transversion rates) of 2 and a strict molecular clock. The Markov chain Monte Carlo (MCMC) was run for 15,000,000 iterations and stored every 1000. Convergence was assessed using Tracer v1.7.2 (Rambaut et al. 2018), where all ESS values were over 500. TreeAnnotator was used to combine iterations, after discarding the first 10% as burn-in, and FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/) was used to visualise the tree.

As there are limited fossil calibration points available within a reasonable phylogenetic distance from Eurasian otters, to provide effective estimations of intraspecific divergence times, priors were selected based on estimated divergence times of our outgroup (L. sumatrana) from our species (L. lutra), taken from previous large scale dated phylogenies. Of the four previous studies to estimate the divergence time of L. sumatrana and L. lutra, two scenarios are supported by two studies each. Firstly, Koepfli et al. (2008) estimated that L. sumatrana and L. lutra diverged at 1.8 million years ago (Mya, 0.9–2.7 95% Highest Probability Density, HPD) based on eight fossil calibration points across all mustelids, including a Lutra affinis fossil from 3.6 Mya (early Pliocene). Secondly, Law et al. (2018) estimated the same species diverged at 1.67 Mya (0.99–2.5 95% HPD) based on 74 fossil calibration points across all Musteloidea species. Thirdly, Waku et al. (2016) estimated much earlier divergence times occurring at 3.24 Mya (1.80–4.81 95% HPD), based on two fossil calibration points >14 Mya, and finally, Hassanin et al. (2020) estimated this divergence to be 3.4 or 3.8 Mya (using either log-normal or uniform fossil calibration priors, respectively) based on 74 fossil calibration points across Carnivora. Based on these studies, we set the time of divergence between L. sumatrana and L. lutra using an external prior, to obtain intraspecific divergence times within L. lutra. Specifically, we used two alternative calibration time scenarios for the L. sumatrana and L. lutra divergence, occurring at (1) 0.9–2.7 Mya and (2) 1.80–4.81 Mya, implemented in BEAST2 as uniform priors, and we present results from both scenarios.

Cytochrome b

Cytochrome b sequences were gathered from GenBank and combined with all whole mitochondrial genome sequences used in this study, totalling 91 sequences from 141 samples. Sequences were aligned and trimmed to 1140 bp of the cytochrome b based on the annotated reference mitochondrial genome (NC_062277.1). A TCS network, and summary statistics were generated using the same methods as for the whole mitochondrial genome data. Full methods for cytochrome b analyses are provided in Supplementary Material. For sample KOR02, uncertainty between an A or G base at 282 bp made it not possible to determine as either haplotype KO3 or KO4, and it is therefore handled as a separate haplotype for summary statistics, rather than incorrectly classifying as either KO3 or KO4, but we acknowledge that it is not a unique haplotype.

Results

Based on the alignment of 27 sequences of 16,432 bp, 26 unique haplotypes across 890 segregating sites were identified. For the cytochrome b region, 141 sequences were included for 1140 bp, from which 27 unique haplotypes were identified. Known cytochrome b haplotypes made up 13 of these, and the remaining 14 haplotypes were novel or unnamed, and therefore we named them as LUT01 to LUT14 (see SM1).

The TCS network of the whole mitochondrial genome sequences identified five main lineages (Figure 1(a)), the geographic locations of which are visualised in Figure 2. Although less clear due to the limited resolution provided by the shorter fragment, the cytochrome b TCS network results (Figure 1(b)) were consistent with the findings from the whole mitochondrial genome sequences, showing the same five main lineages (geographic locations presented in Figure SM 1).

Figure 1. Statistical parsimony networks of Eurasian otters based on (A) whole mitochondrial genomes and (B) cytochrome b. (A) Network of 27 whole mitochondrial genome sequences (16,432 bp, following removal of the tandem repeat region), where nodes are labelled by sample ID/accession and sample location (see SM1). (B) Network of 141 cytochrome b sequences (1140 bp), where nodes are labelled by haplotype name (see SM1). The main five mitochondrial lineages identified in this study are circled and labelled by lineage name, and subspecies found in that lineage. Nodes are a pie chart, colour coded by the number of samples from each sample location, grouped into broad geographic regions, and sized according to haplotype frequency. Dashes on lines represent mutational steps.

Figure 2. Map of sample locations, sample counts, subspecies classifications and whole mitochondrial genome lineages. Exact sample locations are represented by a pin, and where this was not available, the centroid of a region is also represented by a pin (details of exact and centroid locations in SM1). The number of samples from a location is given in white within the pin if it is more than 1. The subspecies classifications of each sample are given in black text beside the pin, with adjacent, duplicate subspecies labels removed. The full species range is highlighted in orange hash from (Roos et al. 2015). Of the 27 Eurasian otter samples included in this study, the locations of 22 are plotted, 1 is unknown, and the only location information given for 2 was ‘Korea’ and for 2 was ‘China’, therefore these are not included in this map.

For the whole mitochondrial genomes data, when samples were grouped by geographic origin, the number of segregating sites and nucleotide diversity were significantly higher in the Japanese samples due to both Lineages 1 and 5 being present (Table 2). Korea showed the lowest nucleotide diversity, followed by China and Western Europe, with sample sizes of 6, 6 and 5. Regardless of sample sizes (n = 2), 2–3× higher nucleotide diversities were observed in Russia and Southeast Asia than Korea, China, and Western Europe. This is likely due to the large geographic range captured by the 2 Russian samples, and the 2 distinct lineages captured by both the Southeast Asian and Russian samples. When samples were grouped by lineage, the lowest nucleotide diversity and number of segregating sites were observed within Lineage 3 (primarily identified in Western Europe) despite its broad geographic range. Lineage 2, found in samples from Korea and east Russia, showed higher nucleotide diversity and segregating sites, and even higher values were found for Lineage 1 (primarily identified in Asia). Due to only one haplotype being identified for both Lineage 4 and Lineage 5, mitochondrial diversity within these lineages could not be determined.

Table 2. Summary statistics of whole mitochondrial genome diversity of Eurasian otters calculated across broad geographic regions and mitochondrial genome lineage. Sample size (n), number of haplotypes (nh), segregating sites (s), haplotype diversity (h) and nucleotide diversity (π). Southeast Asia includes samples from Laos and Sumatra, and Western Europe includes samples from UK, Denmark and Norway.

	N	nh	s	h	π
Samples grouped by broad geographic region:
Korea	6	5	113	0.93	0.0015 (0.0004)
China	6	6	118	1.00	0.0028 (0.0007)
Western Europe	5	5	134	1.00	0.0032 (0.0009)
Japan	2	2	443	1.00	0.0207 (0.0147)
Russia	2	2	116	1.00	0.0070 (0.0050)
Southeast Asia	5	5	208	1.00	0.0052 (0.0014)
Unknown	1	1	–	–	–
Samples grouped by whole mitochondrial lineage:
Lineage 1	13	13	170	1.00	0.0025 (0.0004)
Lineage 2	7	6	113	0.95	0.0015 (0.0003)
Lineage 3	5	5	29	1.00	0.0007 (0.0002)
Lineage 4	1	1	–	–	–
Lineage 5 (L. l. nippon)	1	1	–	–	–
All samples:	27	26	890	1.00	0.0064

Summary statistics for the cytochrome b data reflect the patterns observed in the whole mitochondrial genome sequence summary statistics (see Supplementary Table 1). Despite many fewer samples included for whole mitochondrial genomes (n = 27) compared to cytochrome b sequences (n = 141), a similar number of unique haplotypes were identified; 26 whole mitochondrial genome haplotypes compared to 27 cytochrome b haplotypes. The number of segregating sites (s) and haplotype diversity (h) are lower in the cytochrome b dataset than whole mitochondrial genome sequences, reflecting the differences between the markers. Nucleotide diversity is comparable between datasets, with the exception of a decrease from 0.005 in whole mitochondrial genomes to 0.003 in cytochrome b for Southeast Asia due to the addition of many samples of the same haplotype identified in Laos. A similar effect is observed for China.

For the dated phylogeny, priors and estimated mean divergence times (and 95% HPD intervals) of each lineage identified in Figure 1 are given in Table 3. The resulting topology estimated from the whole mitochondrial genome phylogenetic analysis was consistent between the two tested scenarios (which assumed two different divergence times between L. lutra and L. sumatrana; Figure 3). As expected, dates of divergence among lineages reflected the differences in priors between scenarios, with all divergence dates being earlier when the prior for outgroup divergence was earlier (scenario 2). In line with previous studies, lineage 5 (L. l. nippon sample) was the earliest to diverge within this dataset, between 0.29 and 1.55 Mya, i.e. within the Pleistocene. The 95% HPD intervals overlapped for the remaining divergences among lineages, and all divergences occurred prior to 0.08 Mya.

Figure 3. Dated maximum clade credibility tree from BEAST2 of Eurasian otter whole mitochondrial genome lineages. Topology, branch lengths and 95% HPD interval bars for divergence times were estimated using priors from scenario 1, which assumed a divergence between L. lutra and L sumatrana at 0.9–2.7 Mya. Values at nodes show estimated mean date of lineage divergence based on priors from scenario 1 (left value)/scenario 2 (right value). Time is given in units of Million years ago (Mya). All node posterior probabilities were 1, except the divergence of Lineage 4, which was 0.86 for scenario 1 and 0.87 for scenario 2.

Table 3. Calibration priors and estimated divergence times of Eurasian otter lineages as estimated using BEAST2. All units are Million years ago (Mya), and in brackets are 95% HPD intervals.

	Scenario 1	Scenario 2
Calibration prior
L. lutra/L. sumatrana divergence prior (min–max) in Mya	0.90–2.70	1.80–4.81
Estimated time of divergence in Mya (95% HPD)
L. lutra/L. sumatrana divergence	1.48 (0.90–2.46)	2.83 (1.80–4.43)
Lineage 5 (L. l. nippon) divergence	0.51 (0.29–0.85)	0.97 (0.59–1.55)
Lineage 3 divergence	0.22 (0.12–0.37)	0.42 (0.25–0.67)
Lineage 4 divergence	0.21 (0.12–0.35)	0.40 (0.23–0.64)
Lineage 1 & 2 divergence	0.15 (0.08–0.25)	0.28 (0.16–0.45)

Discussion

Using a combination of newly generated and previously published whole mitochondrial genome and cytochrome b data, we here provide a range-wide synthesis and timeline of the phylogeography of Eurasian otters. While we acknowledge the limited sample size, distribution and subspecies captured in this dataset, we believe this work provides a significant development to our understanding of the species and subspecies relative to prior studies.

Mitochondrial lineages and their discordance with current subspecies

Whole mitochondrial genome Lineages 1 and 2 were originally defined by Waku et al. (2016), and subsequently their geographic distribution was expanded by du Plessis et al. (2023), and further expanded in this study. Specifically, here we identified Lineage 1 to include two additional samples from China (from Hong Kong and Kinmen Island) and three additional samples from Laos. We also identified three additional samples from South Korea in Lineage 2. This confirms the known geographic distributions of these lineages: Lineage 1 is found in China, Japan, Laos and the UK, and Lineage 2 is found in Korea and east Russia. du Plessis et al. (2023) described a previously uncharacterised lineage (Lineage 3), to which the present study assigns two additional samples from the UK, demonstrating the presence of this lineage in western Europe and central Russia.

In addition, we here describe a further two mitochondrial lineages of Eurasian otters: Lineage 4 containing the distinct L. l. barang sample from Sumatra, and Lineage 5 containing the distinct L. l. nippon from Japan as identified in Waku et al. (2016). We found no evidence of lineage 5 in Korea nor neighbouring mainland Asia, however more sampling would be required to investigate if the (extinct) Japanese L. l. nippon is present in Korea (our first hypothesis).

Of the five samples labelled as L. l. barang, one from Sumatra (SRR5192239) and four from Laos (LC094961, OP554563, OP554564, OP554565), the Laos samples were assigned to Lineage 1, and only the sample from Sumatra was identified as Lineage 4. Similarly, the two samples labelled as L. l. nippon (LC049955 and LC050126) were not grouped together but separated into Lineages 1 and 5. Although L. l. chinensis samples are only found in Lineage 1, L. l. lutra samples are found across Lineages 1, 2 and 3. Overall, there is therefore no concordance between the subspecies classification based on Hung and Law (2016) and whole mitochondrial genome lineages, supporting our second hypothesis.

The cytochrome b results highlight the difficulties in identifying the distinct lineages observed from the whole mitochondrial genomes when assessing only a 1140 bp fragment. However, using samples sequenced for cytochrome b increased the sample size and geographic breadth assessed in this study. When comparing summary statistics of the whole mitochondrial genome and cytochrome b datasets, nucleotide diversity was highlighted as a robust metric for use in comparisons among different mitochondrial fragments, showing much less variation between datasets than the number of segregating sites or haplotype diversity. Extracting the cytochrome b region from all whole mitochondrial genome sequences included in our analyses enabled us to identify corresponding lineages across both datasets. Firstly, the sample size of the divergent Lineage 5 (L. l. nippon) increased to two when combining whole mitochondrial genome and cytochrome b data, with the additional sample also being from Japan. Furthermore, we have expanded the likely geographic range of Lineage 3 to include southern Europe (Italy and Portugal), and Iraq. Specifically, three haplotypes (H1, H2 and H3) were identified in samples from Iraq, two of which were also identified in whole mitochondrial genome samples from central Russia, UK and Norway, all from Lineage 3. Samples from Iraq were from two locations, one in northern Iraq, near to the range of L. l. lutra (containing haplotypes from Lineages 1, 2 and 3), and one from southern Iraq, close to or within the distribution range attributed to L. l. meridionalis, for which the whole mitochondrial genome has not yet been sequenced. Our results tentatively indicate that this subspecies (L. l. meridionalis) is likely to belong to Lineage 3, consistent with possible inclusion in L. l. lutra. This result highlights not only the importance of populations in Asia for their disproportionately large genetic and taxonomic variation relative to the well-studied western European populations, but also the importance of populations across the range, for example in the Middle East and North Africa, for which demographic information is difficult to obtain, much less samples for genetic or genomic analysis.

This discordance among subspecies classifications and mitochondrial lineages has been identified in other Eurasian mammals, such as the Eurasian lynx (Lynx lynx) (Mengüllüoğlu et al. 2021; Behzadi et al. 2022), in contrast to species such as the red fox (Vulpes vulpes) which shows considerable overlap between subspecies and mitochondrial lineages in North Africa and the Middle East (Statham et al. 2014). These examples, and the Eurasian otter, highlight the growing need for combined genetic and morphometric analyses to confidently delineate species, subspecies or distinct lineages from intricate species and subspecies complexes, as demonstrated by Colella et al. (2021) for ermines (Mustela spp.).

Diversification timeline and regional genetic diversity

Overall, whole mitochondrial genome haplotype diversity was high; we only observed one duplicate haplotype, highlighting the value of genomics relative to prior genetic methods (Keis et al. 2013). When samples were grouped by location, we observed the highest diversity in regions where more than one lineage was found, for example Japan and Russia. Within lineages, the lowest mitochondrial variation was observed in lineage 3 (found in western Europe), followed by lineage 2 (Korea and east Russia), and the highest variation was observed in lineage 1 (primarily occurring in Asia). These results support our third hypothesis, that mitochondrial variation is higher in Asia (lineage 2 in Korea and lineage 1 in China) than western Europe (lineage 3). These high levels of genetic diversity further highlight the urgent need for genomic studies into the subspecies and populations of Eurasian otters in Asia.

Across Eurasian mammals, high mitochondrial diversity and multiple distinct lineages have been identified in Asia (Li et al. 2021). Most notable are distinct Japanese lineages found in multiple species, for example in Eurasian otters (Waku et al. 2016; Park et al. 2019), brown bears (Ursus arctos) (Hirata et al. 2013), black bears (Ursus thibetanus) (Kishida et al. 2022; Zou et al. 2022), red foxes (Vulpes vulpes) (Statham et al. 2014; Kutschera et al. 2013), and grey wolves (Canis lupus) (Koblmüller et al. 2016). These studies have suggested that multiple waves of recolonisation have resulted in diverse mitochondrial lineages in Japan, with Koblmüller et al. (2016) highlighting the role of sea level changes in the geographic separation of specific lineages. For red foxes, these geographically restricted populations in Japan, show lower genetic diversity than the geographically widespread lineages (Kutschera et al. 2013), in contrast to what we have observed in the Eurasian otter.

The obtained dated phylogeny is consistent with prior studies (Waku et al. 2016), but provides a higher resolution of recent divergences within Eurasian otters than previous work. Divergences among Lineages 1, 2, 3 and 4, are likely to have occurred during the Late Pleistocene, at least 0.08 Mya, whereas the divergence of Lineage 5 is likely to have occurred earlier. However, due to a paucity of suitable fossil calibration points for Eurasian otters and close relatives, it is difficult to provide narrower confidence intervals of divergence times for this species. Assuming these mitochondrial lineages evolved in geographic isolation, separated around 0.08–1.55 Mya depending on the calibration prior scenario, it is likely that regions containing multiple lineages arose through recent natural or artificial translocations. Incidentally, both regions identified in this study to contain multiple mitochondrial lineages in close geographic proximity, Japan and UK, have been hypothesised to have been impacted by recent translocations of the species. du Plessis et al. (2023) identified signals of recent translocations of Asian Eurasian otters to Britain, based on both nuclear and mitochondrial genomic data, and Nakanishi and Izawa (2019) hypothesised that Japan is currently being recolonised by otters from the Korean peninsula, based on camera trapping data. Unfortunately, without modern samples from Japan we cannot validate this hypothesis, however we have included the sample LC049955 published by Waku et al. (2016), which is from Japan, from 1915 or 1916, which groups within Lineage 1 (and mainly comprises samples from China and Laos). Assuming that Lineage 5 (L. l. nippon; represented by sample LC050126) has been present in Japan since at least the 1900s (likely very much longer), the presence of Lineage 1 (sample LC049955) suggests that Lineages 1 and 5 may have coexisted in Japan at least throughout the 1900s, rather than a recent recolonisation of Japan by Lineage 1 postdating the L. l. nippon extinction in 2012. Future work should characterise samples from Japan and surrounding regions over time and to include nuclear data where possible, in order to understand the presence and potential coexistence of multiple mitochondrial lineages in Japan.

This work adds to a growing list of studies investigating range-wide phylogeography of mammals across Eurasia, identifying common and unique biogeographic processes which have shaped the mitochondrial diversity observed in modern populations. Here, we provide the first whole mitochondrial genome assessment of subspecies classifications of Eurasian otters, with a focus on subspecies and populations in Asia. We find discordance between genetic lineages and subspecies, and therefore recommend further research to fully resolve and identify the genetic lineages in the region. We also recommend the use of the identified lineages and nomenclature, to inform conservation and management of the species in Asia, for example when discussing translocations among populations, to avoid anthropogenically mediated gene flow of evolutionarily distinct and isolated populations.

Acknowledgements

We would like to thank all present and past Cardiff University Otter Project staff and volunteers, and the Kadoorie Farm and Botanic Gardens for the collection, storage and donation of a sample (CHN01). The Bioinformatic analyses were performed on the Cardiff School of Biosciences’ Biocomputing Hub HPC/Cloud infrastructure, and we are immensely thankful to Dr Ian Merrick and Andy Ells for their support. We thank the editor and two anonymous reviewers for their helpful feedback. Author contibutions: SJdP, FH, and SH designed research, BL and SH collected samples and coordinated sequencing, SJdP conducted analyses and wrote paper with guidance from FH and K-PK. All co-authors provided feedback on the paper.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Unique haplotype data have been deposited to NCBI Nucleotide Database with the accession numbers OR655422-OR655427 and BK064847. Cardiff University Otter Project has licences to possess and transport otters for England, Wales and Scotland. Permission to possess and transport sample CHN01 was obtained from the Agriculture, Fisheries and Conservation Department of Hong Kong.

Additional information

Funding

This study was supported by the funding of SH by the National Research Foundation of Korea [NRF-2021R1C1C2004162] and the UK Natural Environment Research Council supported this work through the GW4+ Doctoral Training Partnership [award reference: 2194854] of SJdP. The collaboration with K-PK was enabled by a Global Wales International Research Mobility Fund Award [UNIW/RMF-CU/16] to SJdP.