Introduction

Describing and explaining patterns of genetic population structure continues to be a central theme in evolutionary biology (Richardson et al., 2016). Large genetic variation within and among populations is assumed to be important for adaptive evolutionary change, and in this context gene flow among populations may potentially enhance or hamper local adaptation (Wright, 1943; Slatkin, 1987; Hendry et al., 2000). Understanding the origin and evolution of population genetic structure in natural populations is inherently challenging as variation in e.g. population size, life cycle, generation time and behavior all may be important factors impacting on structure (Nadler, 1995; Barrett et al., 2008).

Most studies of population genetic structure tend to focus on single species, but clearly the structure of many organisms depends on how other species around them are distributed and structured. This should be particularly so for parasites, in particular for those that are highly host specific. For instance, ecto-parasites experience a habitat that may be temporary and varying in quality as the host moves across environments, whereas endo-parasites may experience a more stable environment (Poulin, 2007). Contrasting the genetic structure of the host and its ecto-parasite may illuminate how the genetic structure of the parasite is contingent upon its main host (Huyse et al., 2005; Poulin, 2007; Mazè-Guilmo et al., 2016). However, little is known about the specific mechanisms that facilitate or constrain gene flow in parasite taxa (Poulin, 1997; Blouin et al., 1999). Intuitively, host dispersal is important for facilitating gene flow in parasites (Blouin et al., 1995; McCoy et al., 2012). Parasites with only freshwater hosts (termed an autogenic life cycle) have lower dispersal ability between freshwater localities than parasites transmitted through an alternation between freshwater and terrestrial hosts (allogenic life cycle) (Criscione & Blouin, 2004; Blasco-Costa et al., 2013). How the genetic variation and structure covary in associated hosts and parasites seem to be species-specific and dependent upon type of life cycles (Raeymaekers et al., 2008; Louhi et al., 2010; Song et al., 2014). However, remarkably few studies have compared the genetic structure of fish hosts and their parasites (but see Geist & Kuehn, 2008; Sprehn et al., 2015; Kmentová et al., 2020).

Salmonid fishes are commonly hosts to a variety of ecto-parasites (Pickering & Christie, 1980; Pike & Lewis, 1994). Salmonids often have complex life histories that include ontogenetic habitat shifts and changes in social behavior (Hendry et al., 2004; Hesthagen & Sandlund, 2004). Migration between spawning habitats suitable for egg and juvenile development and feeding habitats suitable for larger fish are common (Hesthagen & Sandlund, 2004; Quinn, 2005; Jonsson & Jonsson, 2011). Natal homing is also typical for salmonids (Stabell, 1984), leading to genetic differentiation of neighboring populations even when there are no physical barriers to dispersal (Fleming & Reynolds, 2004; Waples, 2004). Salmonid fishes often change from being solitary and territorial to moving in groups (schooling). In particular, most salmonids will actively aggregate during spawning (Fleming & Reynolds, 2004).

The European grayling (Thymallus thymallus (L.)), hereafter: grayling, is assumed to be a philopatric species, meaning that it returns to its natal spawning location for reproduction (Kristiansen & Døving, 1996). Such natal homing behavior lead to reduced gene flow among populations and thus genetic structuring (Barson et al., 2009; Koskinen et al., 2001; Junge et al., 2011). During early spring, mature grayling living in lakes aggregate outside the spawning tributaries before they enter the tributary for spawning (Beauchamp, 1990; Northcote, 1995). During the spawning season the males defend territories and during such periods of physical aggregation individuals will be close together and parasite transmission may be more common than during the rest of the life cycle (Pickering & Christie, 1980; Richards et al., 2012). If this is the case, and if homing is strong, one may expect the population genetic structure of host and parasite to be associated. The parasite is here highly contingent upon host homing behavior.

Grayling is the host of the ecto-parasitic monogenean Gyrodactylus thymalli Žitňan, 1960 (Soleng et al., 1999). G. thymalli has an autogenic life cycle where transmission occurs directly between grayling individuals when in close contact (see Bakke et al., 2007). Inter-host contacts are most likely occurring during the spawning time of the grayling (Pickering & Christie, 1980). During the remaining time of the year there will probably be less opportunity for physical contact among grayling individuals facilitating parasite transmission. However, it is documented that Arctic grayling (Thymallus arcticus) in rivers establish dominance hierarchies facilitated by individual interactions (Hughes, 1992). We do not know if such hierarchies are developed for grayling feeding in the littoral zone of lakes. Thus, there is most likely a short spatio-temporal parasite transmission window during the spawning season. If this is the case, it will make the Gyrodactylus thymalli-grayling system a good model system for exploring how genetic structure develops among parasites and hosts.

We have compared the population genetic structure of G. thymalli (using mitochondrial DNA sequencing) and its main host, the grayling (using microsatellite genotyping), in the largest Norwegian lake – Lake Mjøsa. Lake Mjøsa holds numerous tributaries where populations of grayling migrate for spawning during spring (Kristiansen & Døving, 1996). Previous studies have detected significant genetic structuring at large geographic scales for both grayling and G. thymalli (Junge et al., 2014; Pettersen et al., 2015; Van Leeuwen et al., 2018). Here, we investigate the genetic structure at a smaller geographic scale, by contrasting samples of hosts and parasites sampled during the grayling spawning season in ten tributaries to Lake Mjøsa. There are no physical barriers to gene flow among tributaries.

Methods

Study area and sampling

Lake Mjøsa is a large (362 km2) and deep (mean depth 153 m, maximum depth 453 m) lake in southeast Norway (L’Abée-Lund et al., 2009) (Fig. 1). The catchment area is 16,563 km2, containing mountainous areas with glaciers as well as lowland areas heavily influenced by human settlement and agriculture (Haugen et al., 2008). The fish assemblage is relatively rich in a Norwegian context, containing 20 species in total (Hesthagen & Sandlund, 2004). Grayling are common throughout the lake and are known to spawn during spring in a large number of tributaries draining into the lake. A more detailed description of the tributaries is given in Kristiansen & Døving (1996). Grayling in Lake Mjøsa usually become sexually mature at an age of three years and mature grayling migrate to the tributaries for spawning during May and June (Kristiansen & Døving 1996). Grayling usually spawn when the water temperature reaches approximately 5°C (own observations and Fabricius and Gustafson 1955; Parkinson et al., 1999; Gregersen et al., 2008). Generally, the grayling complete spawning within approximately two weeks, and then return to the lake. Grayling embryos hatch in June and the juveniles stay in the tributaries until late summer (Kristiansen & Døving, 1996).

Fig. 1
figure 1

Tributaries in Lake Mjøsa where spawning grayling (Thymallus thymallus) and Gyrodactylus thymalli were sampled. The colour codes refer to relative frequencies of the mtDNA (NADH5) haplotypes of G. thymalli (G08–G31). The tributary codes refer to the names in Table 1. In parentheses are given first the total number of sampled grayling and secondly the number of infected grayling (bold). The inset shows a map of Norway with the location of the Lake Mjøsa indicated by a star.

Full size image

In total, 146 mature grayling varying in size from 27 to 44 cm were sampled from ten tributaries by electrofishing during the 2008 and 2009 spawning seasons (Fig. 1, Table 1). Five of the tributaries drain into the eastern side of the lake and four into the western side. The last population was sampled in the main inflowing river (Gudbrandsdalslågen), downstream of the Hunderfossen waterfall (label: HUN) in the north. The geographic distances between pairs of tributaries are available in Additional file 1. There are no apparent physical barriers to migration and gene flow among any of the tributaries. A previous mark-recapture study has shown that grayling tend to spawn repeatedly in the same tributary and that they only to a limited degree migrate between the eastern and western side of the lake (Kristiansen & Døving, 1996). We therefore grouped the tributaries into three different units (west, east and north; Table 1) when testing for large-scale geographic population structure. In addition, we tested for differentiation among the individual spawning populations.

Table 1 Summary table for the number of sampled grayling and Gyrodactylus thymalli, followed by locality code and name and its location within the Lake Mjøsa (west/east/north)
Full size table

All fish were killed quickly by a blow to the head, according to national regulations. All the fins were excised and stored in 96% ethanol for future examination. A small piece of fin tissue from each individual grayling was stored for later DNA-extraction and genotyping. In the lab, the fins of each fish were inspected for the presence of G. thymalli using a stereo-microscope. In total, 129 G. thymalli specimen were sampled for DNA-extraction and sequencing (Table 1). Species of Gyrodactylus are viviparous (Malmberg, 1993) and most or all individuals on one host are likely the offspring of one mother. Thus, we sampled one parasite from each fish, rather than several parasites from the same fish, to maximise sampled genetic variation. A subset (n = 45) of the sampled parasites have been used in a previous study investigating the genetic structure on a much larger geographic scale (Pettersen et al., 2015). All G. thymalli sampled from one tributary are referred to as belonging to one host population.

DNA from fin tissues of the host and whole specimen of the parasite were extracted using the GeneMole® tissue kit (MoleStrips™ DNA) and apparatus (Mole Genetics AS), following the manufacturer’s instructions.

Molecular analyses of G. thymalli

All G. thymalli specimens were sequenced for the mitochondrial dehydrogenase subunit 5 (NADH5) gene. This mtDNA gene has been identified as a highly variable region in G. thymalli (Plaisance et al., 2007; Pettersen et al., 2015). We used the protocol and primer pairs presented in Huyse et al. (2008) to amplify an 894-bp fragment of the NADH5 gene. Polymerase chain reaction (PCR) was performed using PuReTaq ready-to-go PCR beads (GE Healthcare), 1 μmolL−1 of each primer, and 5 μL of the extracted DNA in a 50-μL reaction volume. PCR products were purified using 10 × diluted exoSAP-IT (USB). Cycle sequencing, using the same primers as in the PCR reaction, was performed in 10-μL reactions using 2-μLBigDye terminator cycle sequencing ready reaction kit (Applied Biosystems), 2 μL 5 × sequencing buffer, 10 pmol primer, and 3-μL cleaned PCR product. The cycle sequencing was done on an ABI 3730 high-throughput capillary electrophoresis machine. All sequences were edited using Sequencher 5.0 and aligned with the ClustalW algorithm in MEGA 5.0 (Tamura et al., 2007), with default parameter settings. Genetic distances were calculated according to the Tamura (1992) three-parameter + Gamma (T92 + G) model. This model gave the best fit in a model selection test of 24 different nucleotide substitution models, selected with models with lowest BIC scores (Bayesian information criterion) in MEGA 5.0 (Tamura et al., 2007; Tamura et al., 2004). For each codon, we estimated the number of synonymous and non-synonymous substitutions produced using maximum likelihood (ML) analysis of natural selection codon-by-codon using HyPhy in MEGA 5.0 (Felsenstein, 1981; Muse & Gaut, 1994; Pond & Frost, 2005; Suzuki & Gojobori, 1999). Further, Tajima`s neutrality test was performed on each population in MEGA 5.0 to test if sequences deviated from neutrality expectations (i.e. equilibrium between mutation and genetic drift; (Tajima, 1989). Deviations from neutrality are indicated by p-values below 0.05 in MEGA 5.0.

An analysis of molecular variance (AMOVA) for haplotypes was conducted in GenAlEx 6.2 (Excoffier et al., 1992; Peakall & Smouse, 2006; Peakall & Smouse, 2012) partitioning the variance components among tributaries, among individuals within each tributary, and between tributaries on the eastern and western side of Lake Mjøsa. Additionally, GenAlEx was used for an analysis of pairwise genetic differentiation (population haploid pairwise genetic differentiation φST analog to FST) between tributaries. The φST was then used in a principal component analysis (PCA) to visualize potential G. thymalli groupings (Peakall & Smouse, 2006).

In order to test if there were significant differences in the number of sampled haplotypes among the three regions (east, north and west) a sample-based rarefaction analysis in Estimate S (Version 8.2.0) was used correcting for sample size. To infer the relationships among haplotypes a median-joining network was estimated in NETWORK 4.5.1.6 (Flexus Technology). MEGA 5.0 was used to test for the phylogenetic relationships among haplotypes using Neighbor-Joining (NJ) (Saitou & Nei, 1987) and ML with a bootstrap consensus tree inferred from 1,000 replicates (Felsenstein, 1985) using the same molecular evolution model as before (Tamura three-parameter + Gamma (T92 + G)). Two sequences were used as outgroups; G. thymalli from the type locality (River Hnilec, Slovakia; GenBank number EF527269) and a more phylogenetically distant species G. derjavinoides Malmberg, Collins et al., 2007 (GenBank number EU293891) sampled from brown trout (Salmo trutta L.).

Molecular analyses of grayling

All grayling individuals were genotyped for a set of 17 microsatellite loci following Junge et al. (2010). Methodological details and marker information are found in Additional file 2. Briefly, we ran eight PCRs (using 1 × Qiagen Multiplex PCR Master Mix) with annealing temperatures between 58°C and 60°C. Subsequently, PCRs were combined in two running sets, with three and five PCRs, for electrophoresis on an ABI3730xl Genetic Analyzer, Applied Biosystems (ABI, California, USA). Genotypes were scored using GeneMapper 4.0 software (ABI).

We tested for neutrality of the 17 grayling microsatellites using Lositan (Antao et al., 2008; Beaumont & Nichols, 1996). GenAlEx 6.41(Peakall & Smouse, 2012) was used to estimate grayling microsatellite diversity (frequencies and mean of alleles per locus, observed and expected heterozygosity). To calculate allelic richness the program FSTAT 2.9.3.2 (Goudet, 1995) was used, with the smallest population size set at 8 individuals. Deviation from Hardy–Weinberg and linkage equilibrium was tested using GENEPOP version 4.0.7 (Rousset, 2008). Markov Chain (MC) parameters chosen for the test were: 10,000 dememorization number and maximum number of iterations per batch for 1,000 batches. To avoid type II errors we applied sequential Bonferroni corrections (Moran, 2003).

To determine the genetic structure, a number of different approaches were used. First, an AMOVA was performed using the locus-by-locus procedure in Arlequin 3.5.1.3 (Excoffier et al., 1992; Excoffier et al., 2005) to test for genetic variation across tributaries and/or regions (west/east). Here, the variance components were partitioned into (i) among regions (west/east), (ii) among tributaries and (iii) among individuals within tributaries. The analysis was conducted for the complete dataset, but excluding the northern population HUN, as it was a priori expected to be genetically distinct from the others. The program POWSIM was used to assess the statistical power, by mimicking sampling from populations at a predefined level of expected divergence through random number computer simulations under a classical Wright-Fisher model without migration or mutation. Significance of estimates was based on 1,000 independent simulations.

The program GENEPOP (Weir & Cockerham, 1984) was used to test for genetic differentiation among each population pair at every locus and over all loci by an exact G tests. As described above we used the same MC parameter settings and multiple test correction procedures. We also used GENEPOP to estimate pairwise and global FST (Weir & Cockerham, 1984).

The program STRUCTURE 2.3 (Pritchard et al., 2000) was applied to infer population structuring. Ten replicate runs for K values between 1 and 9, using the admixture model, the LOCPRIOR option (Hubisz et al., 2009) and assuming correlated allele frequencies (Pritchard et al., 2000; Falush et al., 2003) were performed with a burnin of 500,000 and 1,000,000 MCMC iterations. The analyses were performed both for the complete dataset and for the dataset with the HUN sampling location excluded. HUN clearly grouped distinctly from the other populations.

Grayling and G. thymalli

We tested for an isolation-by-distance (IBD) genetic structure (Wright, 1943) in both grayling and G. thymalli using a Mantel test implemented in GenAlEx 6.41 (Peakall & Smouse, 2006; Peakall & Smouse, 2012) by correlating genetic distances (FST/(1 − FST)) or (φST/(1 − φST)) (Rousset, 1997) with geographic distances (km). Two distance metrics were tested, one measure of geographic distances was the shortest water distance between tributary mouths and second measure was the shortest distance along the shoreline (see Additional file 1). These tests were conducted for both grayling and G. thymalli separately. Shoreline distances were used in addition to the shortest distance since mark-recapture studies in grayling in Lake Mjøsa have shown that grayling usually disperse along shorelines rather than across open water (Kristiansen & Døving, 1996). The level of significance was evaluated by performing 9,999 permutations.

In addition, we tested if the pairwise genetic distances between hosts (FST/(1 − FST)) and between parasites (φST/(1 − φST)) were correlated using a Mantel test in GenAlEx 6.41.

Results

G. thymalli

A total of 129 individual G. thymalli were sampled from the grayling spawning in the 10 tributaries to Lake Mjøsa (Table 1). The mtDNA alignment of the NADH5 gene from G. thymalli (768 bp, GenBank no. KY264209-KY264338, Additional file 3) revealed 13 haplotypes, defined by 16 polymorphic sites, leading to 9 amino acid substitutions (Table 1). All the haplotypes have been detected previously in a large-scale study encompassing the total River Glomma watershed (Pettersen et al., 2015). The nucleotide frequencies on NADH5 were 31% (A), 31% (T/U), 16% (C), and 22% (G). The polymorphic sites of NADH5 showed no indication of selection on amino acid triplets based on the codon-by-codon ML analysis.

Sample size, allele frequency (Af), number of segregating sites (S), nucleotide diversity (π), and Tajima`s D test statistics from the different tributaries can be found in Table 1. None of the populations showed divergence from neutrality when tested by Tajima`s D. The two most common haplotypes (G8 and G24) were found in almost all tributaries (Fig. 1). Haplotype diversity was slightly higher on the western side of Lake Mjøsa (8) compared to the eastern side (5). However, the sample-based rarefaction analysis showed overlap in confidence intervals between the number of haplotypes in the western (mean [95% CI]; 10 [6.7–13.3]), eastern (5 [1.7–8.7]) and northern HUN (5 [1.8–8.5]) populations. Private alleles were observed in three different populations; SKU, STA and HUN.

The AMOVA revealed that 94% (VC: 1.233, φST = 0.052, P = 0.030) of the haplotype variance could be explained by within-tributary variation, whereas only 0% (VC: 0.00, φCT = − 0.010, p = 0.680) of the variation could be explained by region (West/East) and 6% (VC: 0.080, φSC = 0.061, P = 0.020) by among-tributary variation.

The network analysis on NADH5 resulted in a star-shaped network with no signal of phylogeographic structuring (Fig. 2). The network analysis detected one unsampled ancestor node (mv1) with four branches with the two most common haplotypes (G24 and G8) in different branches. There was no support for phylogenetic relationships between G. thymalli haplotypes within Lake Mjøsa based on the NJ method and ML method, evidenced by low bootstrap support (i.e. bootstrap estimates below 70%, (Hillis & Bull, 1993).

Fig. 2
figure 2

Median-joining network for the NADH5 haplotypes of Gyrodactylus thymalli. Haplotypes colour-coded as in Fig. 1. The cross lines indicate the mutation steps and the size of the circles represents haplotype frequency. The red node mv1- is the hypothetical unsampled ancestor haplotype

Full size image

We found weak evidence for some pairwise genetic differentiation among G. thymalli sampled in the different tributaries; i.e. 5 out of 45 estimates were significant (Table 2). φST values varied between 0.109–0.764 which is significant given the Bernoulli method. The sample from the GAU tributary was involved in three of these five significant comparisons. There was no relationship between pairwise genetic differentiation (measured as φST/(1 − φST) and pairwise geographic distance measured as either the shortest distance across water (R2 = 0.015, pMantel = 0.636) or the shortest distance along the shoreline (R2 = 0.008, pMantel = 0.220). This means there was no signal of IBD.

Table 2 Distance matrix of population pairwise FST values of grayling (Thymalus thymallus) (below diagonal) and pairwise φST Gyrodactylus thymalli (above diagonal) and result from exact G test (significant pairwise comparisons indicated)
Full size table

A PCA on the population genetic distances revealed that the populations could not be grouped by habitat/geography (Additional file 4).

Grayling

All the microsatellite loci confirmed to neutral expectations when tested and were therefore included in the analyses (Additional file 5). The average number of alleles per locus within a sampling site varied between 4.6 and 6.4 (Additional file 6). Allelic richness (AR) was calculated for all tributaries and ranged from 4.2 to 5.0 (Table 1). Observed (Ho) and expected (He) heterozygosity ranged from 0.56 to 0.64, and 0.60 to 0.67, respectively (Table 1). All populations confirmed to Hardy–Weinberg equilibrium and all locus pairs showed linkage disequilibrium, i.e. they were not linked.

Out of 45 pairwise population comparisons for genetic differentiation, a total of 29 (64%) comparisons showed significant differentiation (Table 2), which is significant given the Bernoulli method (i.e. the probability of getting this result with α = 0.05 by chance alone is very small). The HUN population was significantly differentiated from all the other nine populations. For the remaining samples from Lake Mjøsa tributaries, BAL was significantly differentiated from all other samples except STA (its closest neighboring tributary), whereas FJA, HAM and VIK were differentiated from all but two other tributaries (RIN and STA).

The AMOVA revealed that close to 99% (VC: 5.412, FCT = 0.001, P < 0.5) of the genotype variation could be assigned to within-sample variation, whereas 0% (variance component, VC: 0.058, FST = 0.012, P < 0.01) could be explained by region (west/east) and 1% (VC: 0.064, FSC = 0.012, P < 0.01) by among tributary variation. The global FST (jackknifed over loci) was 0.033 (95% CI 0.004–0.020, 99% CI 0.002–0.024), and pairwise FST-values varied between − 0.008 and 0.094 (Table 2).

There was no clear isolation-by-distance signal, but a weak trend was observed when using the complete dataset and shortest waterway distance as the spatial metric (R2 = 0.530, pMantel = 0.056, Fig. 3; see Additional file 7 for shortest shoreline distances, R2 = 0.014, pMantel = 0.126). After excluding the HUN population (assumed to be a resident river-living population), the relationship between genetic and geographic distance remained non-significant (waterway distance: R2 = 0.003, pMantel = 0.397, shoreline distance: R2 = 0.003, pMantel = 0.335).

Fig. 3
figure 3

Thymallus thymallus: Isolation-by-distance for the complete dataset, based on comparison of pairwise level genetic differentiation with minimum waterway distance between in km (R2 = 0.53, pMantel = 0.056)

Full size image

The STRUCTURE analysis revealed that the most likely number of genetic clusters were K = 2 (Additional file 8), separating the large tributary HUN from the nine smaller tributaries which jointly composed the second genetic cluster (Additional file 9, 2A). Removing the HUN population from the data set and rerunning the STRUCTURE analysis revealed no additional sub-structure within the nine tributaries (Additional file 9, 2B).

Host and parasite comparison

There was no correlation between the pairwise genetic distances of G. thymalli and that of its grayling host (R2 = 0.010, pMantel = 0.470).

Discussion

This study shows a lack of strong genetic structure within both G. thymalli and its host the grayling when comparing samples from ten tributaries to the large Lake Mjøsa in Norway. For grayling, two genetic clusters were evident, one comprising fish from the largest inlet tributary while fish from the remaining nine tributaries formed the second cluster. Weak but significant genetic differentiation among the majority of tributaries was observed for grayling, whereas there was no evidence for structuring in G. thymalli. Further, there was no correlation between the pairwise genetic distances of the hosts and those of the parasites, which might indicate that different mechanisms determine their level of gene flow and thus population structure.

Genetic structure of G. thymalli

A total of 13 mtDNA haplotypes were found in G. thymalli in Lake Mjøsa. Most of these haplotypes were shared among parasites found in the ten grayling spawning populations, except for a few endemic haplotypes. Further, there was no clear geographic structure, neither when tested among the three geographic groups (east–west-north) nor among tributaries. All the haplotypes have been observed earlier in a large-scale study of G. thymalli in the River Glomma watercourse (Pettersen et al., 2015). Lake Mjøsa is part of this river system. The genetic network that was found in Lake Mjøsa is very similar to the network described for the total River Glomma system (Pettersen et al., 2015). Thus, none of the haplotypes in G. thymalli have likely arisen post-glacially in Lake Mjøsa.

There was low support for a geographic structuring of G. thymalli haplotypes among tributaries to Lake Mjøsa based on the Neighbor-Joining method. This was also supported by the EstimateS analysis comparing haplotype diversity adjusting for sample size, and by the AMOVA analysis showing that most of the genetic variation was due to the within-tributary effects. The PCA plot of genetic distances also supported a lack of differentiation. These results all imply that there is a low degree of genetic differentiation among G. thymalli sub-populations in Lake Mjøsa.

No similar studies of genetic variation of G. thymalli in lake ecosystems can be found in the literature. However, a study of 309 G. thymalli from 20 sites in the large Glomma river system found 35 distinct NADH5-haplotypes (Pettersen et al., 2015). Further, Lindqvist et al. (2007) observed three G. thymalli NADH5-haplotypes on a sample of 16 grayling in one site in a tributary. In contrast to the finding of several haplotypes of G. thymalli within different tributaries in our study, only one single haplotype of G. salaris and G. thymalli was found on several populations of Atlantic salmon (Salmo salar L.) and grayling in northern Baltic rivers (Kuusela et al., 2009). Kuusela et al. (2009) related their findings to the fugitive strategy hypothesis (Platt & Weis 1985). This hypothesis suggests that when a river is already occupied by a parasite strain (or haplotype), there is reduced probability for other strains to invade (Kuusela et al., 2009). This situation is related to the concepts of invasion dynamics, the probability of successfully invading an already occupied habitat, and the waiting time until dominance of the invader (Gavrilets & Gavrilets, 2004). There is no support for the fugitive strategy hypothesis in our extensive study. A potential explanation for a one-tributary-one-haplotype observation can be low sampling effort (Walther et al., 1995), as increased sample size leads to an increased number of haplotypes observed (Gavrilets & Gavrilets, 2004).

Genetic structure of grayling

The results suggest that two well-defined genetic clusters of grayling are present in Lake Mjøsa. Grayling sampled in the largest tributary to Lake Mjøsa (HUN) were clearly differentiated from grayling spawning in the nine other tributaries. This finding was also supported by the observed isolation-by-distance trend where the genetic distinctiveness of the HUN tributary was evident. Grayling in the HUN tributary could belong to a stationary river-living population, in contrast to tributary-spawning but lake-living grayling in the other tributaries. This could reduce the probability for gene flow between grayling populations belonging to these two clusters. Furthermore, the HUN tributary is by far the largest in terms of water flow and probably population size, meaning that gene flow from other tributaries likely would have minimal effect on the population structure of grayling in the HUN tributary. However, some gene flow from HUN to other tributaries is expected as indicated in the STRUCTURE plot (see Additional file 7). Junge et al. (2011) studied the population structure of grayling in Lake Lesjaskogsvatnet, located upstream the locality HUN. By studying genetic differentiation among these populations over several years they observed significant among-tributary genetic differentiation even if the level of differentiation was low. This is similar to our observations in Lake Mjøsa. Thus, the grayling in the two lakes are similar in partitioning most genetic variation to within tributaries. Further, there was a weak isolation-by-distance signal when all the populations were included in the analysis. However, the signature became non-significant when the HUN-population was excluded. In comparison, Junge et al. (2011) detected a significant isolation-by-distance pattern for grayling spawning in different tributaries to the much smaller Lake Lesjaskogsvatnet. This could be due to the large environmental differences between tributaries in Lake Lesjaskogsvatnet, potentially leading to stronger selection for natal homing.

The low genetic differentiation seen among the Lake Mjøsa grayling was a little surprising as an earlier mark-recapture study documented moderate to strong homing behavior of grayling in Lake Mjøsa (Kristiansen & Døving, 1996). Further, that study also found that grayling mainly used the littoral zone nearby the tributary where they were initially marked when feeding during summer (Kristiansen & Døving, 1996). The same study found that straying among spawning tributaries was 16% (an overall estimate), and only 2% migration between the western and eastern side of the lake. A recent telemetry study in Lake Lesjaskogvatn has documented a relatively small home range for lake-living grayling (mean 0.58 km2) (Bass et al., 2014). Further, a study of the feeding and habitat use showed that grayling mainly used the littoral zone (Haugen & Rygg, 1996). Thus, given the apparent littoral zone preference of grayling and the moderate to strong homing behavior to spawning tributaries in Lake Lesjasjkogsvatnet, we expected a higher degree of genetic differentiation in the Lake Mjøsa grayling.

Genetic studies have also been conducted on brown trout spawning in several tributaries to Lake Mjøsa (Skaala, 1992; Linløkken et al., 2014). Brown trout from the different spawning tributaries were significantly differentiated, with clear difference between east and west populations (Skaala, 1992). Further, there seemed to be a significant isolation-by-distance pattern (Linløkken et al., 2014). Thus, brown trout seems more genetically differentiated than grayling in spawning tributaries to Lake Mjøsa. Brown trout spawn during autumn, with incubation of eggs and embryos during winter and spring. The importance of being locally adapted may thus be stronger, leading to stronger selection for natal homing than in a spring-spawning fish having rapid development during spring and early summer. Also, for grayling spawning in large rivers, low levels of natal homing and genetic structuring have been found, opposite to what is common for brown trout (Junge et al., 2014; Van Leeuwen et al., 2018).

Grayling—G. thymalli comparisons

There was no correlation between the pairwise genetic distances of grayling and G. thymalli sampled in the different spawning tributaries to Lake Mjøsa. Based on the finding that grayling repeatedly spawn in the same tributaries year after year (Kristiansen & Døving, 1996), we expected that grayling populations should be genetically differentiated and that G. thymalli should follow that same pattern. The estimates of pairwise genetic differentiation among spawning tributaries were significant in 64% of the cases in grayling (FST) and in only 11% in G. thymalli (φST). This pattern resembles the pattern found by Song et al. (2014) who compared an acanthocephalan (Acanthosentis cheni Amin, 2005) parasite and its fish host the Japanese grenadier anchovy (Coili nasus Temminck & Schlegel, 1846). They observed a higher level of gene flow among the parasites than among hosts. This may be a general phenomenon in host-parasite systems (Anton et al., 2007; Levin & Parker, 2013; McCoy et al., 2012). However, in contrast, Criscione & Blouin, (Criscione & Blouin, 2007) studied the trematode parasite (Plagioporus shawi (McIntosh, 1939) Margolis, 1972) and its salmonid hosts (Oncorhynchus spp.) and observed congruent genetic structure and divergence between the host and parasite.

Nieberding et al. (2004) suggested that the genetic variation of a parasite can be used to provide more in-depth information on the phylogeographic history of its host. They suggested that the “magnifying glass” hypothesis, suggesting higher genetic variation and differentiation in the parasite compared to the host, could render the study of the genetic structure of parasites useful for investigating the structure of its host (Mackenzie, 1987; MacKenzie, 2002). In our system, however, the host showed lower levels of gene flow among populations than did the parasite. This led to some genetic structuring in grayling, but less in G. thymalli. These results do not support the “magnifying glass” hypothesis, but are similar the findings in another Gyrodactylus - host system (Huyse et al., 2017). That study compared the genetic variation of G. gondae Huyse, Malmberg & Volckaert 2005 and its sand goby host Pomatoschistus minutus (Pallas, 1770) along the Atlantic coasts of Europe, a much larger area than our study system (Huyse et al., 2017). The main findings were that the parasite populations were significantly differentiated, but the structuring was not consistently higher than in its host populations (Huyse et al., 2017). A recent study of monogeneans (Dactylogyridae) parasitize clupeid fishes in Lake Tanganyika also did not find support for the “magnifying glass” hypothesis (Kmentová et al., 2020). Rather, the parasites did not show any clear geographical population structure. Overall, there seems to be little support for the magnifying glass hypothesis.

It is commonly suggested that the genetic structure of parasites is congruent with that of their hosts, because the host dispersal limits parasites dispersion pattern (e.g. Jarne & Theron, 2001, review by Mazè-Guilmo et al., 2016). However, a meta-study covering diverse taxonomic groups of hosts and parasites, concluded that the genetic differentiation on average was lower in parasites than their hosts (Mazè-Guilmo et al., 2016). This result corresponds well with findings in the present study. Genetic differentiation depends on generation time and population size. G. thymalli has a shorter generation time than grayling (Bakke et al. 2007), and thus stronger genetic differentiation may be expected. However, genetic structure is weaker in parasites with asexual reproduction compared to parasites with sexual reproduction (Mazè-Guilmo et al., 2016). G. thymalli basically reproduce asexually, thus leading to lower overall differentiation rates (Fox et al., 1996; Huyse et al., 2005; Liu et al., 1996).

We assumed that the between-host transmission of G. thymalli mainly occurs when grayling is spawning. In our study we thus sampled grayling only during the spawning season, when grayling concentrate on the spawning areas and interact intensively. However, transmission of parasites among hosts can also occur during time periods outside the spawning season, as seen in grayling when they search for food in the littoral area (Eloranta et al., 2011). In a study by Richards et al. (2012) transmission of Gyrodactylus spp. on guppies were dependent upon the social behavior of the fish where the most parasitized fish had the highest shoaling intensity. Grayling may develop dominance hierarchies, at least in rivers (Hughes, 1992). Such social interactions may lead to parasite transmission. Further, if grayling individuals aggregate while feeding on abundant food resources in the littoral zone transmission may also occur. Thus, transmission of G. thymalli can likely occur also outside of the spawning season, facilitating a more effective transmission of the parasite. Studies on the social behavior of grayling outside of the spawning season is needed.

Conclusions

A moderate to high level of gene flow exists among spawning populations of grayling and among the populations of its ectoparasite G. thymalli. This led to some genetic structuring in grayling, but less so in G. thymalli. This is contrary to our expectation of a moderate to strong association between host-parasite genetic structures and genetic distances. The low genetic divergence in G. thymalli could be due to weak homing behavior in grayling or more social interaction between individuals also outside of the spawning season. The exact mechanisms for the transmission of G. thymalli between grayling host individuals will need to be investigated in the future, taking space and time into consideration.