High levels of inbreeding with spatial and host-associated structure in lice of an endangered freshwater seal.

Host-specialist parasites of endangered large vertebrates are in many cases more endangered than their hosts. In particular, low host population densities and reduced among-host transmission rates are expected to lead to inbreeding within parasite infrapopulations living on single host individuals. Furthermore, spatial population structures of directly-transmitted parasites should be concordant with those of their hosts. Using population genomic approaches, we investigated inbreeding and population structure in a host-specialist seal louse (Echinophthirius horridus) infesting the Saimaa ringed seal (Phoca hispida saimensis), which is endemic to Lake Saimaa in Finland, and is one of the most endangered pinnipeds in the world. We conducted genome resequencing of pairs of lice collected from 18 individual Saimaa ringed seals throughout the Lake Saimaa complex. Our analyses showed high genetic similarity and inbreeding between lice inhabiting the same individual seal host, indicating low among-host transmission rates. Across the lake, genetic differentiation among individual lice was correlated with their geographic distance, and assignment analyses revealed a marked break in the genetic variation of the lice in the middle of the lake, indicating substantial population structure. These findings indicate that movements of Saimaa ringed seals across the main breeding areas of the fragmented Lake Saimaa complex may in fact be more restricted than suggested by previous population-genetic analyses of the seals themselves.

INTRODUCTION

Many large vertebrates, particularly birds of prey and terrestrial and marine mammals, are listed as endangered because of habitat destruction, pollution, overexploitation, direct persecution, or climate change (Courchamp et al., 2018; IUCN, 2021). Species belonging to the charismatic megafauna often act as flagship or umbrella species that attract public attention to the generally dire situation of natural ecosystems across the world (Berti et al., 2020; Thompson & Rog, 2019). In many cases, endangered large vertebrates also constitute important model systems for studying the genetic effects of population bottlenecks and habitat fragmentation (Gousy‐Leblanc et al., 2021; Luo et al., 2019; O'Brien et al., 2017). Conservation‐genetic studies on endangered animals have focused on inbreeding and loss of genetic diversity (Karamanlidis et al., 2021; Rey‐Iglesia et al., 2021), both of which can add to the direct threats imposed by reduced population size, such as Allee effects (Courchamp et al., 1999; Nagel et al., 2021) and sensitivity to environmental and demographic stochasticity (DeWoody et al., 2021; Díez‐del‐Molino et al., 2018; Kyriazis et al., 2021; Lande, 1993; Spielman et al., 2004; Williams et al., 2021).

A fact often overlooked is that charismatic megafaunal species themselves constitute the habitat of other organisms. Large vertebrates host a multitude of ecto‐ and endoparasites, including lice, fleas, nematodes, and cestodes (Pérez et al., 2006; Thompson et al., 2018; Vlasman & Campbell, 2004). In particular, highly host‐specific parasites (i.e., those found on only one host species) may be more endangered than their more obviously threatened hosts (Carlson et al., 2017; Dharmarajan et al., 2021; Dunn et al., 2009; Harris et al., 2014; Pérez et al., 2013; Rózsa & Vas, 2015). While parasitic species are often small and visually unappealing to humans, they still constitute a substantial fraction of global biodiversity and an integral part of healthy ecosystems (Strona, 2015; Thompson et al., 2018), and arguably have their own intrinsic value for ecosystem function and nature conservation (Carlson et al., 2020; Gómez & Nichols, 2013; Stork & Lyal, 1993; Windsor, 1997). Hence, preservation of parasite diversity is important for ensuring normal functioning of both ecosystem‐level processes (Kwak et al., 2020; Milotic et al., 2020) and the immune defences of their hosts (Spencer & Zuk, 2016).

While conservation‐genetic studies have predominantly focused on endangered large vertebrates, genetic investigations of their associated parasites are potentially highly useful for both fundamental and applied research. Two aspects are particularly important:

First, reduced population density of hosts will diminish the chances for movements among host individuals by parasites. Transmission rates will be lowered especially for directly transmitted parasite species, that is, those that require close contact between host individuals for successful transmission. Reduced among‐host transmission probabilities are expected to increase population‐genetic structuring of parasites across host individuals (DiBlasi et al., 2018; Orsini et al., 2013; Sweet & Johnson, 2018). Such effects should also be observed as elevated inbreeding within parasite infrapopulations inhabiting single host individuals, as compared to the whole population (Detwiler & Criscione, 2017).

Second, the spatial genetic structure of parasites can inform us about the population structures, movements, and social networks of their hosts (Gagne et al., 2022; Whiteman & Parker, 2005). Directly transmitted host‐specific parasites will in this respect again be most informative, because their genetic composition will in practice contain a record of past direct interactions among host individuals. From a research perspective, a practical benefit is that parasite genomes are often considerably smaller than those of their vertebrate hosts (de Moya et al., 2021; Kapusta et al., 2017; Zarowiecki & Berriman, 2015). These smaller genomes make approaches that leverage genome sequencing for the collection of population‐genomic data from many individuals more cost effective (Johnson, 2019). The short generation times and faster evolutionary rates in parasites may also mean that differences among subpopulations accumulate faster than in their hosts, potentially allowing analyses of population structuring across finer spatial and shorter temporal scales (Johnson et al., 2014; Martinů et al., 2020; Whiteman & Parker, 2005).

Here, we investigated the levels of inbreeding and genetic differentiation in seal lice (Echinophthirius horridus) living on the endangered Saimaa ringed seal (Pusa hispida saimensis), with respect to both individual host and to geographic space (Figure 1). The Saimaa ringed seal is a postglacial relic subspecies of the ringed seal and is endemic to Lake Saimaa in southern Finland (Figure 2a). The current population of circa 400 individuals is in a slow recovery from a severe bottleneck in the 1980s, when seal numbers were down to less than 150 individuals (Kunnasranta et al., 2021). The postglacial isolation of nearly 10,000 years and the recent severe bottleneck have left their mark in the genetic composition of the Saimaa ringed seal population, which is one of the genetically most uniform pinniped populations on the Earth (Nyman et al., 2014; Palo et al., 2003; Peart et al., 2020; Stoffel et al., 2018). Microsatellite‐based genetic analyses have shown that the fragmented shape of the Lake Saimaa complex (Figure 2a), possibly in connection with the low population size, has led to population‐genetic differentiation across the main breeding areas of the Saimaa ringed seal (Valtonen et al., 2012, 2014).

FIGURE 1

(a) Echinophthirius horridus seal louse male (top) and female (bottom) from Lake Saimaa (for both, ventral view on left and dorsal on right). (b) Seal lice on the muzzle of a dead Saimaa ringed seal; the white arrow shows one of three individuals. (c) Saimaa ringed seal female nursing a weaning‐age pup

FIGURE 2

(a) Map Lake Saimaa, with collection sites of seals and their paired lice shown by coloured dots that are labelled with the seal individual numbers. The main basins of the Lake Saimaa complex are separated by broken lines, with area names indicated on the side. The location of the town Savonlinna at the Kyrönsalmi strait is indicated by a red circle. (b) PCoA ordination plot of seal lice based on their genetic similarity. Lice from the same seal are coloured similarly and connected by lines, and dot colours and shadings correspond to those used in (a). Dot shading indicates the main lake area (see legend). Note that lice from northern Saimaa and Haukivesi tend to be located on the right‐hand side of the ordination, while lice from the two southernmost areas are to the left. (c) Relationship between genetic distance and ln geographic distance between individual seal lice in the full data set. (d) Admixture plot for individual seal lice at K = 2. Section heights within bars show the proportion of ancestry attributed to “northern” (blue) and “southern” (orange) ancestry. Louse individuals are denoted below the plot and ordered from the south to north in the left to right direction, the main lake areas are indicated above the plot, and the locations of the borders between them (see a) are indicated by inverted triangles. The location of the town Savonlinna at the Kyrönsalmi strait is indicated by a red circle above the triangle

Seal lice are in many ways ideal for conservation‐genetic analyses of endangered parasites on endangered hosts. They are obligate, strictly host‐specific parasites that are directly transmitted among host individuals (Leidenberger et al., 2007). Louse genomes in general are small (100–200 Mbp; Allen et al., 2017; Baldwin‐Brown et al., 2021; de Moya et al., 2021), and their generation time is an order of magnitude shorter than those of seals (Kim, 1975; Leonardi et al., 2013; Palo et al., 2003). We estimated levels of genetic diversity and inbreeding, as well as the existence of host‐ and space‐related genetic differentiation, in seal lice on Saimaa ringed seals by sampling pairs of lice from 18 seals across the entire Lake Saimaa complex (Figure 2a). Based on phylogenomic and population‐genomic data sets obtained through genome resequencing of the 36 sampled individuals, we investigated whether lice sampled from the same host individual are on average more closely related than lice on different host individuals, and whether lice show signs of inbreeding on the population and host level. We also investigated whether genetic differences among lice are correlated with their geographic distances, and whether lice show differentiation across the main basins of the Lake Saimaa system. Finally, we contrasted the spatial genetic structures of the lice with results from prior population‐genetic analyses of the Saimaa ringed seal (Valtonen et al., 2012, 2014, 2015).

MATERIALS AND METHODS

Sample collection

Pairs of lice were sampled from 18 seal individuals across Lake Saimaa through 2009–2017 (Figure 2a, Table S1). When more than two lice from a seal individual were available, the specimens for sequencing were selected at random from the infrapopulation sample. The sampling covers all major breeding areas of the Saimaa ringed seal, and the distance between the furthest samples is circa 150 km. We note that our sampling and sequencing design aimed at maximizing the number of seal hosts (i.e., louse infrapopulations) rather than lice per seal in order to optimize possibilities for simultaneously inferring signatures of inbreeding as well as host‐associated and geographic differentiation: (i) host effects on louse population structure and inbreeding can be estimated with two lice per host as long as a large enough number of infrapopulations are sampled, and (ii) testing the effects of geographic subdivision benefits from including as many hosts as possible, rather than from sampling many lice per host. Of the seals, nine were found dead (in the figures and tables below, these hosts are denoted by four‐number codes), and nine were pups briefly captured for radio telemetry studies (below denoted by codes with two letters and two numbers) during long‐term seal monitoring programs of the University of Eastern Finland and the Finnish Forest Management Authority (Metsähallitus). Telemetry studies have been approved by the local environmental authority Centre for Economic Development, Transport and the Environment (permit numbers: ESAELY/433/07.01/2012 and ESA‐2008‐L‐519‐254) and the Animal Experiment Board in Finland (permit numbers: ESAVI/8269/04.10.07/2013 and ESAVI‐2010‐08380/Ym‐23).

In addition to the 36 focal seal lice, two additional specimens were sampled and sequenced (Table S1): individual Echor52, likewise from Lake Saimaa, was sequenced with higher coverage and was used for constructing target gene sequences for mapping reads of the focal lice (see below). The other nonfocal specimen (Echor6) originated from a Ladoga ringed seal (P. h. ladogensis) and was used as an outgroup in phylogenomic analyses of the focal Lake Saimaa lice. All lice were collected into 99.5% ethanol in 2‐ml screw‐cap tubes and stored at −20°C, and each specimen was photographed as a voucher prior to DNA extraction.

DNA extraction and genome sequencing

Whole lice were ground up individually in 1.5 ml Eppendorf tubes, and genomic DNA was isolated using the Qiagen QIAamp DNA Micro Kit (Qiagen). The manufacturer's standard protocol was modified so that specimens were incubated in ATL buffer and proteinase K at 55°C for 48 h instead of the recommended 1–3 h, as well as by substituting buffer AE with buffer EB (elution buffer). This was done to ensure maximal yield (greater than 5 ng) of DNA from the extraction. We quantified each DNA extract with a Qubit 2.0 Fluorometer (Invitrogen).

Libraries for shotgun genomic sequencing were prepared from the extracts with Hyper Library construction kits (Kapa Biosystems). The libraries were quantified by qPCR and sequenced using 150 bp paired‐end reads with an Illumina NovaSeq 6000. These libraries were multiplexed to consume approximately 1/96th of a lane each, producing between 18.4 million and 145 million reads per library (Table S1), representing about 28 to 217× coverage assuming a 100 Mbp genome size. The Echor52 reference sample was multiplexed to consume 1/48th of a lane, producing over 148 million total reads. The FASTQ files from the sequence data were generated and demultiplexed with bcl2fastq version 2.20. All steps of library preparation, sequencing, and data file generation were carried out at the Roy J. Carver Biotechnology Center (University of Illinois, Urbana, IL, USA). Raw reads have been deposited in the NCBI GenBank SRA database (Table S1).

Phylogenomic analyses

To get an overview of relationships among the sampled seal lice, we constructed individual‐level phylogenomic trees based on sequences of 1107 single‐copy protein‐coding orthologue genes. For this, we first used aTRAM (Allen et al., 2015, 2017) to assemble the protein‐coding portions of the focal ortholog genes based on amino acid target sequences from the human louse Pediculus humanus (Johnson et al., 2013). With the aTRAM software, we conducted tblastn searches of 148 M total genomic reads from the seal louse Echor52 library, using the 1107 nuclear ortholog genes from the human louse as blast target references. This software then assembles the resulting blast hits locally in an iterative procedure to produce sequences of the ortholog reference genes for the seal louse. These sequences from the seal louse individual (Echor52) then became new reference sequences for read mapping from other seal lice. We mapped libraries of the focal lice to these assembled target gene sequences from Echor52 using a reference‐mapping pipeline script (https://github.com/adsweet/louse_genomes/) and Bowtie2 (Langmead & Salzberg, 2012). After mapping, we sorted the BAM files and created pileup files using samtools version 1.7 (Li et al., 2009). We used bcftools v.1.7 (Li et al., 2009) to call variants and to convert pileup files to VCF files. Sites with sequence coverage less than 5x or greater than 100x, or with Phred quality scores <28 were filtered (Jiang et al., 2019; McKenna et al., 2010; Sweet et al., 2018) using samtools. From these files, we created consensus sequences for each gene from each individual louse using ambiguity coding for variants.

We aligned nucleotides across all individual lice for each gene separately using pasta version 1.8.2 (Mirarab et al., 2015). Using a custom Python script (Skinner et al., 2020), we removed genes that contained fewer than seven individuals, and then masked sites containing ≥40% gaps using trimal version 1.4 (Capella‐Gutiérrez et al., 2009). After filtering from the 1107‐gene reference set, we were left with 1043 genes (with a total alignment length of 1,379,142 bp) that we used for phylogenomic analyses (Virrueta Herrera et al., 2022). With the aligned data, we performed both a phylogenetic analysis of the concatenated supermatrix and a coalescent analysis of gene trees to produce a species tree. All trees were rooted using the aforementioned louse specimen (Echor 6) collected from Ladoga ringed seal as an outgroup. For the concatenated method, we first ran our gene alignments through RAxML v.8.1.3 (Stamatakis, 2014) and then used the resulting reduced alignment files to create a concatenated matrix. We performed a maximum likelihood (ML) analysis in RAxML, based on a GTR + Γ model of substitution and 100 rapid bootstrap replicates. For the coalescent analysis, we first estimated a tree for each gene alignment in RAxML using a GTR + Γ model, and then summarized the results of the gene‐specific analyses as a coalescent species tree using ASTRAL version 4.10.6 (Mirarab et al., 2014), with quartet‐based local posterior probability support for branches (Sayyari & Mirarab, 2016).

Population‐genomic analyses

We constructed separate population‐genomic data sets for estimation of genetic diversity, inbreeding, and population‐genetic structuring due to seal host individuals and geographic location (Virrueta Herrera et al., 2022). First, we combined the individual VCF files from above into a single VCF file using the merge option in bcftools version 1.7 (Li et al., 2009). We then ran the populations program in STACKS version 2.5 (Rochette et al., 2019) to construct a Genepop‐formatted file containing 2523 SNP sites for use in other population‐genetic analysis programs.

We estimated standard population‐level measures of genetic diversity (number of alleles, observed heterozygosity [H O], heterozygosity within populations [H S], total heterozygosity [H T], and corrected heterozygosity [H'T]) using the Genepop‐formatted file in Genodive version 3.03 (Meirmans, 2020). For the level of individual louse, we calculated the inbreeding coefficient (F) and standardized individual heterozygosity (Coltman et al., 1999) using the –het option in VCFtools version 0.1.15 (Danecek et al., 2011) based on the combined VCF file. The number of sites that could be called as homozygous or heterozygous for individual lice ranged from 2984 to 3066 (mean = 3053.3; Table S1).

To test whether the level of genetic diversity is correlated between lice from the same seals, we used mlRho version 2.9 (Haubold et al., 2010) to calculate sample‐specific mean theta (θ), which is defined as the population mutation rate, or θ = 4N e μ, and which can be used as an indicator of heterozygosity and effective population size (Meyer et al., 2012). For this analysis, we converted pileup files generated from Bowtie2 to profile (.pro) files for each individual louse, and then ran mlRho with maximum distance (M) = 0. These files contained between 1,043,646 and 1,324,364 sites (mean = 1,279,978 sites; Table S1). Finally, we plotted the mean θ of the two lice from each infrapopulation against each other and tested for any correlation between the estimates using reduced major axis regression in the lmodel2 (Legendre, 2018) package in R (R_Core_Team, 2021). We also tested for an effect of lake area and seal host individual on mean θ using GLM ANOVA in IBM SPSS Statistics for Windows version 27.0.1.0, with seal individual nested within lake area in the model.

We inferred the structuring effect of seal host individuals (i.e., infrapopulation structure) by estimating genetic self‐similarity and similarity among individual seal lice based on within‐ and between‐individual kinship coefficients (Loiselle et al., 1995) in Genodive version 3.03 (Meirmans, 2020). As a second estimate of differentiation among infrapopulations, we calculated overall F ST in a data set partitioned by seal host individual in Genodive.

We visualized overall genetic similarities among individual lice by principal coordinates analysis (PCoA) in adegenet (Jombart, 2008; Jombart & Ahmed, 2011) in R. The PCoA method seeks the best approximation in reduced space of a matrix of Euclidean distances. Its principal components optimize the representation of the squared pairwise distances between individuals (Jombart, 2016). We then assessed population structure by estimating the ancestry of individual lice using ADMIXTURE version 1.3 (Alexander et al., 2009). We ran ADMIXTURE for K (number of ancestral populations) = 1–10 with the cross‐validation method to test for the optimal value of K. More optimal values of K will show lower cross‐validation error relative to less optimal values.

To investigate spatial genetic differentiation in the seal louse population within Lake Saimaa, we used two methods:

First, we correlated genetic distances among louse individuals to their geographic distances. Because lice from the same seal cannot be considered independent replicates in an isolation‐by‐distance (IBD) analysis, we added the sampling‐site coordinates of each individual into the Genepop file, but then split the file into 10 separate data sets containing only one randomly selected louse per seal individual. The existence of IBD was then tested for each data set in GenePop version 4.7.5 (Rousset, 2008), based on genetic distances estimated based on the â statistic (Rousset, 2000) and ln geographic distances estimated based on the sampling‐site coordinates. Statistical significance of the regression slopes was inferred on the basis of 95% confidence intervals obtained through ABC bootstrapping (Leblois et al., 2003) and Mantel tests based on 10,000 permutations of individual locations. We also performed a corresponding analysis including all 36 lice, but with the minimum geographic distance among individuals set to 10−5, so that lice from the same seal were not included in the estimation of the regression coefficient.

The second test for spatial effects was done with a hierarchical locus‐by‐locus AMOVA performed in Arlequin version 3.5.2.2. (Excoffier & Lischer, 2010). Prior to the analysis, we converted the Genepop‐formatted data file to Arlequin format using the Widgetcon 1.0.0. website (Aydın et al., 2019) and manual editing. In the analysis, we divided the lice into three main areas (Northern Saimaa + Haukivesi, Pihlajavesi, and Southern Saimaa) defined based on the main basins and breeding areas of ringed seals within Lake Saimaa (Figure 2a). This division scheme is slightly simplified from the one used in the analyses of spatial genetic differentiation in Saimaa ringed seals by Valtonen et al. (2012, 2014), because lice from only a single seal from Northern Saimaa were obtained, so we collapsed this sample of two lice into those from the adjacent Haukivesi population. The AMOVA was then performed with infrapopulations (seal host individuals) nested within lake area and including the level of louse individual in the analysis. Statistical significance of the effect of lake area and infrapopulation was determined by 10,000 permutations of seals (infrapopulations) among lake areas and lice among seals within areas.

RESULTS

Phylogenomic trees

The ML phylogeny based on the concatenated alignment revealed a few clear cases in which lice from the same seal host individual were each other's closest relatives (Figure S1a). The tree also showed some indication of lice from the same area being clustered close to each other, but bootstrap support values for groupings were generally very low across the tree, although this is not unexpected given these are individuals of the same species. Interestingly, by contrast, the coalescent ASTRAL tree revealed a clear structuring effect of seal host individual, with 14 out of 18 sampled louse infrapopulations coming out as monophyletic (i.e., the two individual lice from the same seal host individual were each other's closest relatives; Figure S1b). The general pattern of structuring by lake area was likewise more evident in the coalescent tree, although the support for the backbone structure of the phylogeny was weaker than for the clades formed by infrapopulations, which were in many cases strongly supported (Figure S1b). While the ASTRAL coalescent approach has not often been applied to analyses within a single species, this analysis differs from the concatenated approach in that ASTRAL summarizes the collection of individual gene trees in a coalescent framework. In contrast, the concatenated approach finds the most likely tree for all data combined. Individual gene trees may be more likely to follow patterns of inheritance and gene genealogies within a species, which could explain why the tree from the coalescent approach more closely reflects infrapopulation structure.

Overall observed heterozygosity in the focal seal louse population was 0.199 (s.d. 0.004), expected (total) heterozygosity 0.234 (s.d. 0.004), and corrected expected heterozygosity 0.238 (s.d. 0.004). Heterozygosity within infrapopulations was 0.164 (s.d. 0.003). On the level of individual lice, standardized heterozygosity ranged from 0.173 to 0.284 (Table S1).

As expected, higher individual heterozygosity estimates corresponded to lower inbreeding coefficients (0.290 to −0.167; Table S1). Estimates of individual mean θ ranged between 4.68 × 10−4 and 7.99 × 10−4 (Table S1) and were statistically significantly positively correlated between lice from the same infrapopulation (Figure 3; reduced major axis regression r = .556; p = .016). A statistically significant effect of infrapopulation (nested within region) on mean θ was also revealed in the GLM ANOVA (df = 15, F = 2.929, p = .016), but estimates did not differ across the three regions of the lake (df = 2, F = 2.046, p = .164).

FIGURE 3

Correlation between sample‐specific estimates of mean θ (an indicator of heterozygosity and effective population size) of lice collected from the same seal individuals (i.e., same infrapopulation). Dot colours correspond to those used in Figure 2a, labels indicate the seal individual from which the lice were collected, and dot shading shows the lake area (see legend). The red line represents the correlation from a reduced major axis regression, and grey lines represent the confidence limits of the slope

Between‐individual kinship coefficients ranged between −0.205 and 0.478 and were generally highest between lice from the same infrapopulation (Table S2). Self‐similarities ranged between 0.463 and 0.709, with a mean of 0.585 (s.d. 0.062).

The overall F ST among infrapopulations was 0.312, which was statistically highly significantly different from 0 (p < .001). The population structuring arising from the seal host individual is seen also in the PCoA ordination plot, in which lice from the same infrapopulation tended to cluster together (Figure 2b).

The ADMIXTURE cross‐validation analysis returned an optimal K value of 2 with a CV error value of 0.0012 (Figure S2). The analysis at the optimal K = 2 revealed a sharp change in ancestry proportions roughly in the middle of the lake, corresponding to the limit between Northern Saimaa + Haukivesi and the two southern parts of the lake (Figure 2d).

Plotting the genetic distances among lice against their ln geographic distances revealed a classic IBD pattern (Figure 2c). In the statistical analyses of the relationship using ten subsampled data sets, the mean intercept was 0.274 (range 0.178–0.323) and the mean slope parameter 0.045 (range 0.019–0.053). The bootstrapped 95% CI of the slope parameter did not include zero in any of the subsampled data sets (overall lower/upper range 0.012–0.063). p‐values estimated by Mantel tests were highly significant or significant at p = .0002–.009 in nine cases and marginally significant at p = .061 in one case. For the analysis based on the complete data set of 36 lice, the relationship between genetic and ln geographic distance was estimated as â = 0.224 + 0.043x, with the 95% CI of the slope parameter being 0.035–0.050 and p < .0001 in the Mantel test. Due to the shape of the lake (Figure 2a), calculation of direct interindividual distances from sampling‐site coordinates will slightly underestimate the longest distances across the main basins of the lake, but, given the strength of the IBD pattern (Figure 2c), this is unlikely to have an effect on the general relationship.

The aforementioned patterns were largely summarized by the results of the hierarchical locus‐by‐locus AMOVA, which revealed statistically significant differentiation among the three lake areas as well as among infrapopulations within the areas (Table 1). The differentiation among lice within infrapopulations was strongly and statistically significantly negative, which is a further indication of inbreeding within populations of lice from the same seal individual (Table 1).

TABLE 1

Results of the hierarchical locus‐by‐locus AMOVA when individual lice are grouped according to three main lake areas (Figure 2a) and host seal individuals (infrapopulations) within the areas

Source of variation	df	Sum of squares	Variance components	Percentage of variation
Among lake areas	(2)	2653.232	40.78219	12.99
Among infrapopulations within lake areas	(15)	6452.108	67.27029	21.43
Between louse individuals within infrapopulations	(18)	2890.500	−44.45757	−14.16
Within louse individuals	(36)	9002.000	250.26769	79.74
Total		20,997.840	313.86260

Note: The effect of all explanatory variables is significantly different from 0 at p < .0001.

DISCUSSION

Population‐genetic investigations of endangered parasites can inform us about their population size, genetic diversity, and level of inbreeding, all of which have the potential to influence the likelihood of extinction through deterministic or stochastic processes (DeWoody et al., 2021; Kyriazis et al., 2021; Spielman et al., 2004). Importantly, detailed genetic surveys focusing on parasites of endangered vertebrates have wider applied value, as parasite‐specific analyses can also illuminate biological features of their hosts and, thereby, aid in designing actions for conserving both the parasites and their hosts (Gagne et al., 2022; Sweet et al., 2020; Whiteman & Parker, 2005). Here, we used genome‐level data of seal lice living on the landlocked Saimaa ringed seal to gain insights into the population structure of the lice as well as their lake‐endemic hosts. With a population of barely over 400 individuals, the Saimaa ringed seal is one of the most endangered pinniped populations on the Earth (Kunnasranta et al., 2021). Previous studies have shown that genetic diversity of the Saimaa ringed seal population is extremely low in comparison to other seal species (Martinez‐Bakker et al., 2013; Nyman et al., 2014;Peart et al., 2020; Stoffel et al., 2018). In addition, the main breeding areas of the Lake Saimaa complex harbour partially isolated subpopulations (Valtonen et al., 2012, 2014, 2015). Given that the population density of seals within the lake is low and that seal lice require close contact between host individuals for transmission, we expected that lice inhabiting the same seal would tend to be closely related as a result of within‐host inbreeding. Furthermore, we predicted that the low diversity and distinct spatial genetic structure found in the Saimaa ringed seals would be reflected in the genetic composition of their lice. Our phylogenomic and population‐genomic analyses based on whole‐genome resequencing data from 36 lice sampled from 18 seals across Lake Saimaa indeed supported all of these predictions. Below, we discuss our main results and their implications for the conservation of Saimaa ringed seals, their host‐specific lice, and endangered host–parasite systems in general.

Genetic diversity, differentiation among infrapopulations, and inbreeding

The Saimaa ringed seal is genetically highly uniform in comparison to its sister subspecies in the Baltic Sea (P. h. botnica), Lake Ladoga, and the Arctic Ocean (P. h. hispida; Martinez‐Bakker et al., 2013; Nyman et al., 2014; Palo et al., 2003). The loss of diversity is apparently a consequence of a small founding population and long postglacial isolation (Nyman et al., 2014), as well as the severe anthropogenic 20th‐century bottleneck (Peart et al., 2020; Stoffel et al., 2018). In our focal seal lice, overall heterozygosity estimates (H O = 0.199, H'T = 0.238) were not extremely low. In a study by DiBlasi et al. (2018), mean heterozygosities were 0.449 for pigeon body lice and 0.557 for wing lice. However, these latter values were based on microsatellite markers, which tend to have many alleles per locus and, hence, result in high heterozygosity estimates (Sunde et al., 2020). Our individual‐level estimates of genomic diversity are, however, directly comparable to those of Leonardi et al. (2019), who used the same genomic markers to estimate θ values for five species of seal lice infesting Antarctic and Australian seals having very large population sizes. In their study, species‐specific θ estimates based on individual lice ranged between 0.00107 and 0.00367, which is substantially higher than our individual‐level estimates for E. horridus lice within Lake Saimaa (mean = 0.00062). The highest θ estimate in our data set (0.00080) is also lower than the lowest values (range = 0.00087–0.00863) found by Sweet and Johnson (2018) for seven species of chewing lice on New World ground‐doves. The level of genetic diversity of seal lice within Lake Saimaa therefore seems to directly reflect the low population size and genetic uniformity of their endangered hosts.

From the perspective of parasites of large vertebrates, each host individual constitutes a distinct resource “island” (Itescu, 2019; Koop et al., 2014). If the frequency of among‐host dispersal is low in relation to the generation time of the parasites, parasite populations on different host individuals (infrapopulations) will over time tend to become genetically differentiated from each other (DiBlasi et al., 2018; Huyse et al., 2005). Indeed, our phylogenomic trees, between‐individual kinship coefficients, and estimates of among‐infrapopulation F ST's consistently showed that lice collected from the same Saimaa ringed seal individual are on average genetically more similar than are individuals collected at random from the host population. Notably, the population‐genetic differentiation found across lice collected from different seal individuals is not the only aspect that is affected by the fact that seal lice are distributed into distinct infrapopulations: Using the same specimens that were analysed in this study, Doña et al. (2021) found that infrapopulation identity explained a major proportion of the variation in microbiome composition within individual lice.

The structuring imposed by infrapopulations is clearly visualized in the PCoA ordination, in which lice originating from the same seal are generally located close to each other (Figure 2b). According to the hierarchical AMOVA controlling for within‐lake spatial structure, variation among infrapopulations accounts for 21% of the genomic variation in the louse population (Table 1). Our estimated overall F ST among infrapopulations (0.312) is high in comparison to studies on among‐host differentiation in human body and head lice (F ST = 0.048 in both; Leo et al., 2005), pigeon body (F ST = 0.225) and wing (F ST = 0.075) lice (DiBlasi et al., 2018), and feather lice on Galapagos hawks (pairwise F ST = 0.145–0.183; Koop et al., 2014). Unfortunately, we can only make general comparisons among these different louse–host systems because the spatial scale of different studies varies considerably, and the high heterozygosity of microsatellite markers used in previous studies will in theory suppress estimates of among‐population differentiation (Alcala & Rosenberg, 2019; Jakobsson et al., 2013; Meirmans & Hedrick, 2011). Comparative studies have, however, indicated that microsatellite and SNP markers produce roughly similar estimates of population differentiation (Lemopoulos et al., 2019; Sunde et al., 2020). Hopefully, genomic approaches and the gene ortholog SNP based markers applied here will in the future allow more direct comparisons of genetic variation and differentiation measured from different study systems.

Based on previous studies of louse infrapopulations and the biology of seal louse transmission, it was not unexpected to find genetic differentiation among infrapopulations of the Saimaa seal lice. Despite their long coevolutionary history with aquatic mammals, seal lice are still essentially terrestrial organisms (Leidenberger et al., 2007; Leonardi et al., 2013). Therefore, transmission of lice requires direct contact between seals while they are not submerged in water (Kim, 1975). Within Lake Saimaa, lice are probably transmitted mainly between mothers and pups during nursing (Figure 1c) as is the case in other species of seal lice (cf. Kim, 1975; Leonardi et al., 2013). However, close seal‐to‐seal encounters also occur during the early‐summer moulting period, when two or more seals can share resting sites on large lakeside rocks (Biard et al., 2022). Recent observations also indicate that multiple seals can co‐inhabit the same resting lairs that the seals dig into lakeside snowdrifts during winter (M. Kunnasranta, personal observation), so these may provide additional opportunities for louse transmission.

In addition to the inbreeding caused by transmission dynamics, louse infrapopulations on single seals are also presumably quite small, further increasing the level of inbreeding. In a sample of 49 seals in the collections of the University of Eastern Finland, the number of collected lice ranged from one to 32. Seal lice are difficult to collect exhaustively and immature individuals may go unnoticed, but it seems reasonable to assume that infrapopulation sizes range in the tens rather than in the hundreds. Our sample of seal lice indeed showed clear genomic signs of inbreeding. Individual F values are on average slightly positive and the hierarchical AMOVA showed slightly negative estimates for differentiation between individuals from the same seal. Furthermore, the mean of pairwise Loiselle's kinship coefficients within infrapopulations (0.31) exceeds the expected value between parents and offspring or between siblings (0.25), and the mean of self‐similarity (0.58) likewise exceeds the expectation (0.50) in a randomly‐mating population. Interestingly, the level of genomic diversity and inbreeding varies among infrapopulations, because estimates of θ, which is proportional to the effective population size (Haubold et al., 2010), was found to be positively correlated between lice collected from the same seal individual (Figure 3). This variation evidently reflects substantial differences in infrapopulation size and age, but potentially also stochastic immigration of unrelated individuals into small and generally closed louse infrapopulations.

Spatial differentiation

Spatial population‐genetic differentiation in host‐specific parasites is expected to be influenced by the dispersal patterns of their host species, but spatial structuring can be either weaker or stronger than in the hosts (Cole & Viney, 2019; Dharmarajan et al., 2016; Mazé‐Guilmo et al., 2016; McCoy et al., 2005; Sweet et al., 2020). Weaker differentiation is expected if the parasite species also utilizes intermediate hosts or other host species, has a large effective population size in relation to its host, or if it has a complex life cycle with a highly dispersive life stage (Blasco‐Costa & Poulin, 2013; DiBlasi et al., 2018; Solórzano‐García et al., 2021). By contrast, relatively stronger differentiation is the norm if the parasite is host‐specific, directly transmitted, occurs at low prevalences, and has a comparatively short generation time and high mutation rate (Mazé‐Guilmo et al., 2016).

Despite its large size, Lake Saimaa is in fact a labyrinthine watercourse system formed by several main basins connected by narrow straits (Figure 2a). The fragmented structure of the lake has left its imprint in the genetic composition of the Saimaa ringed seal population, which exhibits an isolation‐by‐distance pattern and differences in the frequencies of mitochondrial haplotypes and nuclear microsatellite alleles across the main breeding areas (Valtonen et al., 2012, 2014, 2015). Similar to the patterns found in the seal hosts, our genome resequencing data of lice revealed a parallel isolation‐by‐distance gradient (Figure 2c) and spatial differentiation (Figure 2b) within Lake Saimaa. As a result, lice from all three main areas of our analysis tended to be grouped together in the PCoA ordination (Figure 2b). Inspection of eigenvalues of the ordination axes additionally shows that most of the variation is explained by Axis 1, which largely corresponds to sampling locations in the north–south direction across the lake.

Importantly, our Admixture results reveal that the main division within the focal seal louse population occurs in the middle of the lake, around the Kyrönsalmi strait (Figure 2a,d). Both shores of the strait are currently covered by the town of Savonlinna, with over 30,000 inhabitants. However, the area has had a substantial human population at least since the foundation of the medieval St. Olaf's Castle on an island in the middle of the strait in 1475 (Taavitsainen, 2005). Given that Saimaa ringed seals were actively hunted until their protection in 1955, the growing human population may have essentially stalled seal—and seal louse—migration between the northern and southern halves of the lake for some five to six hundred years. The low signature of northern genomic ancestry in five lice from Southern Saimaa (Figure 2d) might conceivably result from the experimental translocation of a female seal (Phs152) from Haukivesi to the southern parts of the lake in 1992. This move may have led to inadvertent north–south translocation of lice (and, hence, northern genetic variation), as seal Phs152 is known to have reproduced in its new home range, and it was still alive in 2020 (Kunnasranta et al., 2021).

It is noteworthy that the differentiation in seal lice (Figure 2d) appears to be stronger than that estimated for their seal hosts on the basis of mtDNA and microsatellite data by Valtonen et al. (2012, 2014, 2015). The seal population exhibits statistically significant lake‐wide differences in the frequencies mtDNA haplotypes and microsatellite alleles, but microsatellite‐based assignment analyses by Valtonen et al. (2014) produced spatially restricted clusters only if sampling‐site coordinates were used as background data (priors) in the analyses. In addition, the clusters were not strictly area‐specific, so that individuals belonging to most clusters could be found in several areas of the lake. The stronger spatial signal in lice is most probably due to our much larger genome‐level data set, but also to the fact that seal lice can produce several generations per year (Kim, 1975; Leonardi et al., 2013), while the generation time of ringed seals has been estimated at circa 11 years (Palo et al., 2001). Hence, the seal louse population will accumulate spatial genetic differences substantially faster than their seal hosts.

CONCLUSIONS

Our phylogenomic and population‐genomic analyses of host‐specific ectoparasitic E. horridus seal lice from the lake‐endemic and endangered Saimaa ringed seals show that the louse population consists of genetically distinct infrapopulations that differ among seal individuals and experience high levels of inbreeding. Furthermore, comparisons to genome‐level studies from other louse groups suggest that overall genetic diversity within the focal seal louse population is low—a result that seems to parallel the genetic uniformity of the Saimaa ringed seal population (Nyman et al., 2014; Palo et al., 2003). Comparative studies involving seal lice from the Baltic Sea and Lake Ladoga will be needed for inferring the taxonomic status of the E. horridus population isolated within Lake Saimaa, but our results indicate that the population may be genetically at least as distinct as the Saimaa ringed seal, which has evolved into a separate subspecies after becoming landlocked after the last Ice Age (Kunnasranta et al., 2021; Nyman et al., 2014). Further studies are also required for inferring the ecological and evolutionary relevance of reduced genetic diversity in the focal seal lice. While inbreeding and low genetic variation can suppress viability and reproductive success at the level of both individuals (Blomqvist et al., 2010; Kardos et al., 2016) and populations (Ekroth et al., 2019; Spielman et al., 2004), many parasites are known to experience regular cycles of inbreeding due to their biological characteristics (Appelgren et al., 2018; Detwiler & Criscione, 2017; Van Den Broeck et al., 2014). Hence, parasites may be tolerant to the negative effects of inbreeding (Price, 1980), possibly through purging of deleterious genetic variation (Benesh et al., 2014). Inbred hosts have been shown to be more susceptible to parasitism in many species (Cassinello et al., 2001; Coltman et al., 1999; Hoffman et al., 2014), but far less is known about the effects of inbreeding on parasite performance (Forsman, 2014; see also Benesh et al., 2014; Fredericksen et al., 2021). The endemic Saimaa ringed seals and their specialist lice therefore constitute a promising model system for investigating host susceptibility and parasite infectivity in a “coevolutionary cold spot” in which interactions are highly specialized but in which both hosts and parasites have reduced genetic diversity.

Our population‐genomic analyses revealed a distinct genetic discontinuity in the louse population at the Kyrönsalmi strait, which separates the northern and southern halves of the Lake Saimaa complex. Importantly, this division in the seal louse population suggests that the Saimaa ringed seals of the northern and southern parts of Lake Saimaa are more isolated from each other than mtDNA‐ and microsatellite‐based analyses of the seals themselves have indicated. According to our data, the genetic effects may simply not yet have manifested in the seals due to their longer generation time. To make the comparisons between seals and their lice more comparable, the investigations based on mtDNA and microsatellites by Valtonen et al. (2012, 2014, 2015) should be followed up by genome‐level analyses of the seal population in order to obtain a clear view of their spatial differentiation within the Lake Saimaa complex. Overall, our results highlight how genome‐level analyses of parasites can provide a tractable, cost‐effective, and sensitive early‐warning system for detecting host population fragmentation before the genetic effects are evident in their vertebrate hosts.

AUTHOR CONTRIBUTIONS

Kevin P. Johnson, Tommi Nyman and Stephany Virrueta Herrera conceived the study. Mervi Kunnasranta and Eeva Ylinen obtained samples. Stephany Virrueta Herrera and Kevin P. Johnson collected the data. Stephany Virrueta Herrera, Andrew D. Sweet, and Tommi Nyman analysed the data. Tommi Nyman, Mervi Kunnasranta, Eeva Ylinen, and Kevin P. Johnson obtained financial support for the project. Tommi Nyman and Stephany Virrueta Herrera wrote the manuscript, and all authors contributed to editing the manuscript.

CONFLICT OF INTEREST

The authors declare no competing interests.

BENEFIT‐SHARING STATEMENT

Benefits from this research accrue from presenting information on the biology of an endemic and endangered host–parasite system and the sharing of our data and results on public databases as described above. This study complies with laws governing handling of endangered animals (see Section 2).

Figure S1 Phylogenomic trees for the analysed seal lice.

Figure S2. Cross‐validation error in relation to the number of ancestral populations (K).

Click here for additional data file.

Table S1 Collection information for the sequenced seal lice and their seal host individuals, with information on individual‐level heterozygosity, inbreeding, and data (numbers of sites and reads, and accession numbers).

Table S2. Pairwise kinship coefficients between the analysed seal lice.

Click here for additional data file.