INTRODUCTION
The study of genetic diversity and gene flow can provide fundamental insights into the demographic and evolutionary history of populations (Wright, Reference Wright1943; Slatkin, Reference Slatkin1987). In order for a population to adapt to environmental change, genetic variation is required (Slatkin, Reference Slatkin1987). Gene flow due to individuals dispersing among sub-populations may increase genetic variation. However, the opportunity for individuals to disperse depends on the geographical topography and presence of potential barriers to movement. A lack of gene flow may lead to loss of genetic diversity and increased genetic drift (Grenfell et al. Reference Grenfell, Pybus, Gog, Wood, Daly, Mumford and Holmes2004). In general, there are a large number of studies on genetic diversity in a wide range of study systems. However, studies on the genetic variation and geographical distribution of parasites on larger scales are still rare. Host specific monoxenous fish parasites are interesting for studies of genetic diversity at local and regional scales as their occurrence and genetic diversity depend on the distribution and migration of their hosts as well as transmission probabilities. In rivers and streams, these parasites can disperse downstream by passive drift of detached parasites or fish host migration, whereas upstream dispersal depends on host migration (Criscione and Blouin, Reference Criscione and Blouin2004). Such upstream migration may be restricted by barriers (waterfalls, dams) and may over time lead to a reduced genetic variation upstream due to genetic drift and bottleneck events in small populations. Further, the dominating downstream gene flow can over time result in increased genetic diversity in downstream areas.
Fish ectoparasites of the genus Gyrodactylus have a short generation time, give birth and have no specialized transmission stages. In the recent years, several molecular studies has been done on the taxonomy, systematics, phylogeography and genetic variation of species of Gyrodactylus based on mitochondrial DNA sequences (see e.g. Hansen et al. Reference Hansen, Bachmann and Bakke2003, Reference Hansen, Bakke and Bachmann2007a , Reference Hansen, Bakke and Bachmann b ), but the genetic varation within river systems has not been studied in detail. One of the species that has been studied in some detail is Gyrodactylus thymalli Žitňan, 1960 (Monogenea) and several phylogenetic lineages and haplotypes of this species have been found on its main host grayling Thymallus thymallus (L.) in European rivers (Hansen et al. Reference Hansen, Bakke and Bachmann2007a , Reference Hansen, Bakke and Bachmann b ; Lindqvist et al. Reference Lindqvist, Plaisance, Bakke and Bachmann2007; Anttila et al. Reference Anttila, Romakkaniemi, Kuusela and Koski2008; Kuusela et al. Reference Kuusela, Holopainen, Meinila, Anttila, Koski, Zietara, Veselov, Primmer and Lumme2009). Gyrodactylus thymalli occurrs frequently on grayling in the Glomma river system, the largest watercourse in Norway (Mo et al. Reference Mo, Appleby and Sterud1998; Hansen et al. Reference Hansen, Bachmann and Bakke2003, Reference Hansen, Martinsen, Bakke and Bachmann2006, Reference Hansen, Bakke and Bachmann2007a , Reference Hansen, Bakke and Bachmann b ). Hansen et al. (Reference Hansen, Bachmann and Bakke2003, Reference Hansen, Martinsen, Bakke and Bachmann2006, Reference Hansen, Bakke and Bachmann2007a , Reference Hansen, Bakke and Bachmann b ) found 8 haplotypes on a small number of examined grayling from a restricted part of the river.
The topography of the river Glomma is strongly influenced by the land rise after the last glacial period in Northern Europe that ended about 10 000 years ago. Following de-glaciation, the landscape has been influenced by many geological processes such as land rise, river capture and sea-level rise. As a result of these post-glacial geological processes several waterfalls act as barriers to fish upstream migration. This has influenced the immigration history of the various fish species found in the river Glomma (Huitfeldt-Kaas, Reference Huitfeldt-Kaas1918; Andersen and Borns, Reference Andersen and Borns1994; Nesbø et al. Reference Nesbø, Fossheim, Vøllestad and Jakobsen1999; Koskinen et al. Reference Koskinen, Ranta, Piironen, Veselov, Titov, Haugen, Nilsson, Carlstein and Primmer2000; Bernatchez, Reference Bernatchez2001; Østbye et al. Reference Østbye, Bernatchez, Næsje, Himberg and Hindar2005). Koskinen et al. (Reference Koskinen, Ranta, Piironen, Veselov, Titov, Haugen, Nilsson, Carlstein and Primmer2000) suggested that two immigration routes at different times during the post-glacial process have resulted in two different evolutionary lineages of grayling in the Glomma river system. This could also have affected the genetic variation and genetic population structure of G. thymalli in the river system.
In a recent paper Fromm et al. (Reference Fromm, Burow, Hahn and Bachmann2014) suggested that G. thymalli is a junior synonym to Gyrodactylus salaris Malmberg, 1957 described from Atlantic salmon Salmo salar L. However, Fromm et al. (Reference Fromm, Burow, Hahn and Bachmann2014) did not include specimens from the type localities for the two species which should be used in species synonymiation (zoological nomenclature). Furthermore, G. thymalli is not able to survive on Atlantic salmon while G. salaris can survive on grayling (Sterud et al. Reference Sterud, Mo, Collins and Cunningham2002). Thus, we have chosen to use the name G. thymalli for the parasites in this study.
Here, we present the first in-depth analysis of the genetic variation of a naturally occurring Gyrodactylus species, G. thymalli, from the largest river system in Norway. Given the large geographic scale (see below) and complex immigration history of the host species, we here investigate the level of genetic diversity and how it is distributed across this large geographic scale.
MATERIALS AND METHODS
The study system
The river Glomma covers has a large catchment area in Norway (Fig. 1), draining 13% (41 971 km2) of the country. The river system has two major branches – the east branch (upper Glomma) and the west branch (the river Gudbrandsdalslågen). The two branches merge and form the lower main stem of the river (lower Glomma) about 136 km before draining into the sea. The upper Gudbrandsdalslågen branch includes Mjøsa (MJO), the largest lake in Norway (L'Abée-Lund et al. Reference L'Abée-Lund, Eie, Faugli, Haugland, Hvidsten, Jensen, Melvold, Pettersen, Petterson and Saltveit2009).
Fig. 1. The Glomma river system, with the locations where Gyrodactylus thymalli was sampled (the population codes refer to the names in Table 1). The colour codes refer to the concatenated mtDNA genes NADH5(G) and COI(C) haplotypes (see online Supplementary Table 2 and 3). P.H* indicates private haplotypes.
Table 1. The number of sampled grayling (N f ) and number of infected fish with Gyrodactylus thymalli (N i+G.t.). Locality name and code for each population, its location within the river system, and distance from the outlet to the sea (km) is given (see Fig. 1)
Grayling is found throughout the Glomma river system (Hesthagen and Sandlund, Reference Hesthagen and Sandlund2004). Individual migrations over distances up to 154 km has been observed for grayling in this river (Heggenes et al. Reference Heggenes, Qvenild, Stamford and Taylor2006), but upstream migrations are often prevented by natural waterfalls and dams (Østdahl et al. Reference Østdahl, Skurdal, Kaltenborn and Sandlund2002; Heggenes et al. Reference Heggenes, Qvenild, Stamford and Taylor2006). These barriers have led to population differentiation (Heggenes et al. Reference Heggenes, Qvenild, Stamford and Taylor2006; Junge et al. Reference Junge, Museth, Hindar, Kraabøl and Vøllestad2014). Grayling has a life history with ontogenetic habitat shifts, with spawning and early juvenile life occurring in rivers and streams during spring and summer. The juveniles are territorial, but larger and older individuals may aggregate in schools (Northcote, Reference Northcote1995).
Sampling
Grayling were sampled at 20 localities using various fishing methods during 2007–2011 (Table 1, Fig. 1). The sampling sites were selected to cover the distribution of the grayling in the Glomma river system. The average length between collection sites is 42 km, and the difference between the northern population and the southern populations of the watercourse is 506 km.
After capture, fish were rapidly euthanized by a blow to the head. Individual fish <15 cm were preserved in 96% ethanol and examined whole. From larger fish, only the fins were excised, preserved and examined. In the laboratory, whole fish or fins were examined for the presence of G. thymalli under a stereomicroscope. Although many fish carried several Gyrodactylus specimens, only one specimen was sampled from each fish. The rationale behind this was that species of Gyrodactylus are viviparous (Malmberg, Reference Malmberg1993) and most or all individuals on a fish are likely the offspring of one mother. Thus, sampling one parasite from each fish, rather than several parasites from the same fish, was an attempt to avoid pseudoreplication. In this study, all G. thymalli specimens from several fishes sampled at one locality will be referred to as one population.
Molecular methods
DNA was extracted from Gyrodactylus specimens using the Gene Mole® Robot with MoleStrips™ DNA Tissue kit, following the manufacturer's instructions. The mitochondrial dehydrogenase subunit 5 (NADH5) gene and cytochrome oxidase I (COI), the bar-coding gene (Hebert et al. Reference Hebert, Ratnasingham and deWaard2003), were amplified. The protocol and primer pair described by Huyse et al. (Reference Huyse, Buchmann and Littlewood2008) was used to amplify an 894-bp fragment of the NADH5 gene. For the COI gene the protocol and primer pair described by Meinilä et al. (Reference Mäinila, Kuusela, Ziętara and Lumme2002) was used to amplify 820-bp. Polymerase chain reaction (PCR) was performed using PuReTaq ready-to-go PCR beads (GE Healthcare), 1 µmol L−1 of each primer, and 5 µL of the extracted DNA in a 25 µ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 µL BigDye terminator cycle sequencing ready reaction kit (Applied Biosystems), 2 µL 5× sequencing buffer, 10 pmol primer, and 3 µL cleaned PCR product. The sequencing was done on an ABI 3730 high-throughput capillary electrophoresis machine.
Data analysis
All sequences were proof-read and edited using Sequencher 5.0 and aligned using the Crustal W algorithm in MEGA 5.0 (Tamura et al. Reference Tamura, Dudley, Nei and Kumar2007), with default parameter settings. Genetic distances were calculated separately for COI and NADH5 according to the Kimura three-parameter estimates plus Gamma (T92 + G). This model gave the best fit in a model selection test of 24 different nucleotide substitution models, selected based on the corrected Akaike information criterion (Hurvich and Tsai, Reference Hurvich and Tsai1989) in MEGA 5.0 (Tamura et al. Reference Tamura, Nei and Kumar2004). To test for positive selection, we estimated the numbers of synonymous (s) and nonsynonymous (n) substitutions (dN/ds) using maximum likelihood (ML) analysis of natural selection codon-by-codon by HyPhy in MEGA 5.0 (Felsenstein, Reference Felsenstein1981; Muse and Gaut, Reference Muse and Gaut1994; Suzuki and Gojobori, Reference Suzuki and Gojobori1999; Pond and Frost, Reference Pond and Frost2005; Pond et al. Reference Pond, Frost and Muse2005). A neutrality test (Tajima`s D) was used on each population (Tajima, Reference Tajima1989) to test for sequence deviation from equilibrium between mutation and genetic drift. To test for hierarchical population structure, as well as for a geographical pattern of population subdivision, analysis of molecular variance (MANOVA) in the program GenAlex 6.2 was used (Excoffier et al. Reference Excoffier, Smouse and Quattro1992; Peakall and Smouse, Reference Peakall and Smouse2012). Genetic variance (based on φST) of G. thymalli was partitioned among individuals within populations and among populations (φSC) and between regions (φCT, region 1: Elverum (ELV), Holset (HOL), Rena (REN), Åsta (AAS), Osensjøen (OSE), Koppang (KOP), Storsjøen (STO), Alvdal (ALV), Røros (ROR), region 2: Eidsvoll (EID), Mjøsa (MJO), Hunderfossen (HUN), Trettenstrykene (TRE), Vinstra (VIN), Otta (OTT), Vågå (VAA), Dovre (DOV), region 3: Sarpsborg (SAR), Nittedalselva (NIT), Kongsvinger (KON)) using 9999 permutations. Further, the same program was used to calculate pairwise genetic distances (based on φSC) among populations and further used in a principal coordinate analysis (PCoA). The two axes explaining most of the variation in genetic distance were then extracted. PCoA plots give an opportunity to visualize potential population grouping in relation to geography (Peakall and Smouse, Reference Peakall and Smouse2012).
Under migration-drift equilibrium, pairs of populations are expected to exhibit a significant correlation between their genetic and geographic distances, termed ‘isolation by distance’ (IBD, Wright, Reference Wright1943). This means that populations in close proximity to each other should be genetically less differentiated, due to ongoing gene flow between them, compared with populations that are separated geographically. The occurrence of IBD was tested by correlating pairwise genetic distance with hydrographic distance (km) using a Mantel test implemented in GenAlEx 6.41 (Peakall and Smouse, Reference Peakall and Smouse2012). The level of significance was estimated by performing 9999 permutations. To calculate the hydrographic distance between the sampling locations, the shortest distance over water was used (see online Supplementary Table S1).
A median-joining network analysis (Bandelt et al. Reference Bandelt, Forster and Rohl1999) with maximum parsimony calculation (Polzin and Daneshmand, Reference Polzin and Daneshmand2003) of mtDNA haplotypes was performed using NETWORK 4.5.1.6 (http:www.fluxus-engineerin.com) Flexus Technology software. The phylogenetic relationship between haplotypes based on the mtDNA sequences was inferred using the Neighbour-Joining method (Saitou and Nei, Reference Saitou and Nei1987), ML (Felsenstein, Reference Felsenstein1981) and Maximum Parsimony (Felsenstein, Reference Felsenstein1978) analyses implemented in MEGA 5.0 with bootstrap estimates inferred from 1000 replicates (Felsenstein, Reference Felsenstein1985). Two relevant sequences extracted from Gen-Bank were used to compare with our G. thymalli data (G. thymalli (Gen-Bank acc.no EF 527269), G. salaris (DQ988931)), and a sequence from G. derjavinoides was used as an outgroup (EU293891).
To estimate the time of colonization of G. thymalli in the Glomma river system we used data from Hansen et al. (Reference Hansen, Bachmann and Bakke2003, Reference Hansen, Bakke and Bachmann2007a , Reference Hansen, Bakke and Bachmann b ) indicating that G. thymalli from the river Glomma are highly divergent from G. thymalli in a neighbouring river system (Trysilelva). Hansen et al. (Reference Hansen, Bachmann and Bakke2003, Reference Hansen, Bakke and Bachmann2007a , Reference Hansen, Bakke and Bachmann b ) calculated this divergence to be 2·32% based on the Kimura 2-parameter model (Kimura, Reference Kimura1980). These two lines have most likely not had contact after the post-glaciation immigration of grayling. The time of separation between these two lineages was estimated using the function Inferring Ancestral Sequences ML in MEGA 5.0. NADH5-sequences were not available from the river Trysilelva, therefore only COI sequences were used in this analysis. To be able to compare with previous result (Meinilä et al. Reference Mäinila, Kuusela, Ziętara and Lumme2004; Hansen et al. Reference Hansen, Bakke and Bachmann2007a ), the Kimura 2-parameter model was applied in this analysis and the following sequences from GenBank were used: G. thymalli (AY146612, AY146613 and AY486544) from Trysilelva, and Gyrodactylus lavareti (AY225306) as outgroup. This analysis is based on precise estimates of a molecular clock (i.e. mutation rates) (Kumar, Reference Kumar2005). We used three different mutation rates, one conservative estimate of 2·0% divergence per million years generally suggested for mtDNA (Irwin et al. Reference Irwin, Kocher and Wilson1991; Bermingham et al. Reference Bermingham, McCafferty and Martin1997; Bernatchez, Reference Bernatchez2001), and one intermediate mutation rate of 13·7% and one high mutation rate of 20·3% divergence per million years suggested by Meinilä et al. (Reference Mäinila, Kuusela, Ziętara and Lumme2004) based on studies of G. salaris and G. thymalli.
RESULTS
In total 687 grayling were examined for the presence of G. thymalli; 309 G. thymalli individuals were collected and analysed (Table 1). Gyrodactylus thymalli was found throughout the Glomma river system, with prevalence varying from 16 to 100% (Table 1).
The molecular analysis was based on a final alignment of 777 and 770 bp for the COI and NADH5 genes, respectively. The NADH5-sequences exhibited in total 35 distinct haplotypes (online Supplementary Table S2), defined by 38 polymorphic sites. These mutations resulted in 15 amino acid substitutions (dN/dS = 0·39). The COI-sequences exhibited 37 distinct haplotypes, defined by 28 polymorphic sites, leading to 3 amino acid substitutions (dN/dS = 0·07, online Supplementary Table S3). The nucleotide frequencies in NADH5 were 30·5% (A), 31·3% (T/U), 22·4% (C), and 15·9% (G) and on COI are 28·4% (A), 28·4% (T/U), 21·6% (C) and 21·6% (G). The nucleotide diversity (π) for the two genes ranged between 0·001 and 0·004 in the populations (online Supplementary material; Table S2 for NADH 5 and Table 3S for COI). The ML analysis of natural selection codon-by-codon was not significant for any of the polymorphic sites. Further, there was no divergence from neutrality tested by Tajima's D. For the concatenated NADH5 and COI sequences a total of 58 haplotypes were found, and 80% of the populations had private alleles (Fig. 1, Supplementary Table S2 and S3). The MANOVA using the concatenated NADH5 and COI sequences indicated that 39% of the variation was found among regions, 28% among populations and 33% within populations (MANOVA P = 0·01, Table 2). A principal coordinate analysis (PCoA) of population genetic differentiation (φST,) and the Nei genetic distance showed a clear geographic structure with three main clusters (Fig. 2). One cluster consisted of all populations in the west branch (OTT, VIN, DOV, TRE) plus the ROR sample in the northern part of River Glomma. The second cluster contained populations found in KON and in the lower Glomma (river NIT, SAR near the outlet). The largest cluster contained all the populations from the east branch and from Lake MJO (Fig. 2).
Fig. 2. Plot of the principal coordinate axis 1 (58%) and 2 (28%) of pairwise population genetic distances of Gyrodactylus thymalli populations within the Glomma river system. The analysis is based on the concatenated mtDNA genes NADH5 and COI. The population codes as in Table 1.
Table 2. Analysis of molecular variance (AMOVA) is based on concatenated mtDNA genes NADH5 and COI of Gyrodactylus thymalli. Among regions; are tested under random permutation of individuals across regions (region 1: ELV, HOL, REN, AAS, OSE, KOP, STO, ALV, ROR, region 2: EID, MJO, HUN, TRE, VIN, OTT, VAA, DOV, region 3: SAR, NIT, KON, the population codes refer to the names in Table 1). Among populations; are tested under random permutation of individuals across populations. Within population; are tested under random permutation of individuals across populations without regard to either of their original populations
There was a highly significant correlation between hydrographic distance, and population genetic differentiation (φST, Mantel test, r = 0·28, P = 0·01, Fig. 3), with the exception of a few outliers exhibiting high waterway distances and low genetic differentiation (φST). The ROR population had haplotypes for the NADH5 gene that cluster to the west branch rather than the downstream Glomma river (see Fig. 2). A re-analysis, excluding the ROR population, gave a stronger correlation between waterway distance and φST (r = 0·48, P = 0·01).
Fig. 3. Correlation between hydrographic distance (waterway distance, km) and Gyrodactylus thymalli population pairs genetic distances (φST/1 − φST) within the Glomma river system (Mantel test, r = 0·28, P = 0·01). The analysis is based on the concatenated mtDNA genes NADH5 and COI. The outliers from the Røros population down to the right. A re-analysis, excluding the Røros population, gave a stronger correlation (r = 0·48, P = 0·01).
The geographic distribution showed more or less the same structure for both NADH5 (G), and COI(C). In the west branch the haplotype G14C16 was the most common, compared with the east branch, where the haplotypes G24C20 and G24C29 dominated. Haplotype G04C36 dominated in the lower Glomma (NIT, SAR) and at KON (Fig. 1).
The network analysis detected 11 possible unsampled ancestral nodes. All the three haplotype clusters were linked to one unsampled ancestral node (mv10) in three different branches (Fig. 4). This network indicates a clear geographic structure. In the lower Glomma; haplotype G04C36 grouped with G08C09. From G04C36 a star-shaped network structure grouped the less common haplotypes. All the haplotypes from west branch was found with a star-shaped network structure around haplotype G14C16. The rest of the haplotypes had more or less a star-shaped structure linked with G24C20 and G24C29. There was weak support (bootstrap estimates <80%) for phylogenetic relationships between the haplotypes within Glomma river system when Neighbour-Joining method, ML and Maximum Parsimony analyses were used.
Fig. 4. Median-joining network is based on the concatenated mtDNA genes NADH5(G) and COI(C) haplotypes of Gyrodactylus thymalli, with haplotypes colour coded as in Fig. 1 (For the private haplotypes codes see online Supplementary Tables 2 and 3). The cross-lines indicate mutation steps more than one and the size of the circles represent haplotype frequency. The nodes mv1-11 is hypothetical unsampled haplotypes.
Based on analysis of the COI-gene, the time of separation of the G. thymalli populations in the river Trysil and those in the river Glomma was estimated to be 5000 ± 1000 (95% CI) years ago using the conservative mutation rate (2·0% divergence per million years), 2200 ± 600 years ago when using the intermediate mutation rate and 680 ± 80 years ago using the high mutation rate suggested by Meinilä et al. (Reference Mäinila, Kuusela, Ziętara and Lumme2004) for the genus Gyrodactylus (13·7–20·3% divergence per million years).
DISCUSSION
The monogenean ectoparasite G. thymalli and its grayling host have a wide distribution in the Glomma river system. Infected grayling were sampled at 20 localities and based on the analysis of two mitochondrial genes (COI and NADH5) we found large genetic variation with a high number of G. thymalli haplotypes throughout the river system. There were large genetic differences between populations and these differences were correlated with hydrographic distance. The 37 COI haplotypes revealed in our study included all the 8 haplotypes described earlier (Hansen et al. Reference Hansen, Bachmann and Bakke2003, Reference Hansen, Martinsen, Bakke and Bachmann2006, Reference Hansen, Bakke and Bachmann2007a , Reference Hansen, Bakke and Bachmann b ). By using the concatenated NADH5 and COI sequences, 58 haplotypes were found. The three most common haplotypes were differentially geographically distributed in the river system and represented nodes in a star-shaped network linked with several less-common haplotypes. These three main clusters of haplotypes were linked to three regions (the east branch, the west branch and lower Glomma) and described 39% of the molecular variation. This genetic structure may have evolved from one ancestral population after the last glaciation, or it may be the result of colonization by genetically different G. thymalli populations during separate invasion events. It may be impossible to resolve these two explanations, as earlier phylogenetic analyses were unable to resolve the relationship between the various clusters in Europe (Hansen et al. Reference Hansen, Martinsen, Bakke and Bachmann2006, Reference Hansen, Bakke and Bachmann2007a , Reference Hansen, Bakke and Bachmann b ; Meinilä, et al. Reference Mäinila, Kuusela, Ziętara and Lumme2004). The estimated time of colonization of G. thymalli into the Glomma river system was after the last glacial period in Northern Europe. The different estimates for the divergence of the Glomma G. thymalli from the Trysil G. thymalli varied from 680 to 5000 years ago, varying from very recent divergence to several thousand years after the last glacial maximum.
Current available genetic methods have uncovered that the last ice age has strongly affected the genetic population structure of many species in the Northern hemisphere (Hewitt, Reference Hewitt2000). A high haplotype diversity and low nucleotide diversity, together with a star-shaped haplotype network (single step mutations from central haplotypes), are indicative of a classic post-glacial expansion (Hewitt, Reference Hewitt1996). This fits with our observations for G. thymalli in the Glomma river system. However, the distribution can also be due to different waves of post-glacial immigrant coming from different glacial refugia. Several studies propose that Mid-Norway and Sweden is such a mixture zones for many species (Hewitt, Reference Hewitt2000). Several evolutionary lineages of grayling have been identified, and these lineages seem to have invaded Norway during different time period and following different migration routes (Koskinen et al. Reference Koskinen, Ranta, Piironen, Veselov, Titov, Haugen, Nilsson, Carlstein and Primmer2000; Gum et al. Reference Gum, Gross and Kuehn2005, Reference Gum, Gross and Geist2009). It has been suggested that at least two grayling lineages meet in Lake MJO, the largest water body in the Glomma river system. This is supported by the observation of high numbers of G. thymalli haplotypes, and that all the three haplotype clusters are represented. The grayling and associated G. thymalli can have dispersed upstream or downstream, depending on how water flow has varied through time.
Dispersal and gene flow between populations of host specific fish parasites with a direct life cycle depends on dispersal of the fish host (Criscione and Blouin, Reference Criscione and Blouin2004). Barriers to fish upstream migration will reduce the two-directional genetic exchange between parasite populations and result in separated populations (Blasco-Costa et al. Reference Blasco-Costa, Waters and Poulin2012). The Glomma river system has numerous waterfalls and dams acting as migration barriers for grayling. When grayling colonized the river early after the last glaciation the geographic configuration of the river was different than today. However, in general upstream migration was prevented in a ‘stepwise’ manner as waterfalls arose when the ice melted and the land rose due to isostatic rebound (Huitfeldt-Kaas, Reference Huitfeldt-Kaas1918; Andersen and Borns, Reference Andersen and Borns1994).
For the entire river system, there was a strong correlation between geographic and genetic distance suggestive of a classic isolation-by-distance genetic structure, indicating restricted gene flow between the G. thymalli populations. In this large river system gene flow is limited by the presence of numerous migration barriers such as natural waterfalls and man-made dams. Such restricted gene flow might cause mutations to be sustained in isolated populations that have been separated for a long time (Mills, Reference Mills2007). Our results show that 80% of the populations carry private haplotypes and in two of the populations examined we only observed such private haplotypes (AAS, OSE). These two populations are probably small, and are located where only downstream gene flow seems possible. Small populations of G. thymalli might lead to a strong genetic drift and an increased potential for allele fixation through genetic drift. Low genetic diversity might also be due to inbreeding (Mills, Reference Mills2007). Isolated populations with private haplotypes as found in Glomma have also been found in a study of the genetic structure of the digenean Crassicutis cichlasomae (Razo-Mendivil et al. Reference Razo-Mendivil, Vazquez-Dominguez and de Leon2013). Their findings were that multiple colonization events and subsequent isolation are likely factors that shaped the genetic structure of the parasite (Razo-Mendivil et al. Reference Razo-Mendivil, Vazquez-Dominguez and de Leon2013). The grayling and the associated G. thymalli most likely colonized the mid parts of the Glomma river system early after the ice melted following the last glacial period. Following deglaciation when also the upper parts became available, the land also rose due to the isostatic rebound. At the same time the sea level rose, but at a different rate. Only grayling invading the river early were then able to colonize the upper reaches of the river system. These upper reaches later became isolated, allowing different mtDNA lineages of G. thymalli to evolve in isolation. Later, when the lower parts of the river, including the lake MJO, emerge from the sea and became accessible for fresh water fish, both upstream and downstream migrating grayling could colonize these areas (Koskinen et al. Reference Koskinen, Ranta, Piironen, Veselov, Titov, Haugen, Nilsson, Carlstein and Primmer2000). Koskinen et al. (Reference Koskinen, Ranta, Piironen, Veselov, Titov, Haugen, Nilsson, Carlstein and Primmer2000) found three distinct lineages of grayling mtDNA in Scandinavia, presumably originating from three glacial refugia. Two of these lineages are found in MJO, indicating that grayling from two different evolutionary lineages live in sympatry in MJO. This can probably explain the occurrence of the most frequent G. thymalli haplotypes from all the three genetic clusters in the lake MJO. In addition, several G. thymalli haplotypes were found in the lake only. This observation can be the result of the relatively high sample size that we have from this lake. Or these private haplotypes may have evolved in geographically separated grayling populations spawning in the numerous tributaries draining into MJO (Kristiansen and Døving, Reference Kristiansen and Døving1996).
The genetic diversity of G. thymalli that we observed can be considered as high, whereas the genetic diversity of the grayling host is rather low (Koskinen et al. Reference Koskinen, Ranta, Piironen, Veselov, Titov, Haugen, Nilsson, Carlstein and Primmer2000). Earlier studies have shown similarity in genetic diversity of a fish host and its parasites (see e.g. Wu et al. Reference Wu, Wang, Xi, Xiong, Liu and Nie2009; Atkinson and Bartholomew, Reference Atkinson and Bartholomew2010; Razo-Mendivil et al. Reference Razo-Mendivil, Vazquez-Dominguez and de Leon2013) and parasites can be used as a proxy for understanding the evolutionary history of the host (Nieberding and Olivieri, Reference Nieberding and Olivieri2007). In addition to the natural dispersal of fish and parasites following deglaciation, anthropochore relocation of grayling, within and/or between river systems, may have resulted in dispersal of G. thymalli haplotypes. The historically confirmed anthropochore relocation of grayling from OTT river to Lake VAA in 1906 (Huitfeldt-Kaas, Reference Huitfeldt-Kaas1918) is reflected in the genetic similarity of the G. thymalli haplotypes in Lake VAA and those in the OTT river as well as the genetic similarity between grayling above and below the migration barrier in the system (Junge et al. Reference Junge, Museth, Hindar, Kraabøl and Vøllestad2014). Human relocation of fish might also explain the genetic similarity between the haplotypes G03C16 and G14C16 in the east branch (ROR, ALV) and west branch (VAA, DOV, OTT, VIN, TRE), respectively. This was also supported by the isolation-by-distance analysis where ROR was removed, which showed a stronger overall relationship than when ROR was included. In a population genetic study of grayling from the east branch of Glomma, the ROR population was more closely related to grayling from the neighbouring river Trysil than to populations downstream in the Glomma (ELV, ALV) (Heggenes et al. Reference Heggenes, Qvenild, Stamford and Taylor2006). This result can be explained by a transient bypass between Glomma and Lake Femunden in the river Trysil in 1762. Fish movement to the Glomma from Lake Femunden is documented (Berg, Reference Berg1986).
There is an ongoing discussion about if mtDNA is a good neutral marker, as selection on mtDNA genes has been demonstrated (Ballard and Pichaud, Reference Ballard and Pichaud2014). Our results show that mutations in the NADH5 gene cause amino acid substitutions to a greater extent than shown in previous studies of G. thymalli (Huyse et al. Reference Huyse, Buchmann and Littlewood2008). The differences on the NADH5 gene have often been associated with selection. However, a test for selection (using the ML model) was not significant. The COI gene contained less information at the population level than the NADH5 gene, even if both genes have similar estimated mutation rates. The NADH5 gene has also earlier been reported to have more variation than the COI gene in Gyrodactylus (Lindqvist et al. Reference Lindqvist, Plaisance, Bakke and Bachmann2007; Huyse et al. Reference Huyse, Buchmann and Littlewood2008), and has been identified as being among the most variable regions of the mitochondrial DNA (Plaisance et al. Reference Plaisance, Huyse, Littlewood, Bakke and Bachmann2007). By using both genes, we achieved stronger signals of population structure and a distinct isolation-by-distance genetic structure. In future studies, it is important to increase sample sizes, at the same time as using more sensitive molecular markers for more detailed studies of population structure. Since the whole genome of G. salaris is mapped (Hahn et al. Reference Hahn, Fromm and Bachmann2014), these genetic resources clearly provide more opportunities for detailed studies into population structure by using markers such as single nucleotide polymorphism and microsatellites.
In conclusion, the observed mtDNA gene diversity of the G. thymalli populations in the Glomma river system is the result of postglacial processes that have created barriers to upstream migration of fish leading to isolation by distance structure. Reduced upstream gene flow, as well as potentially high levels of genetic drift, has resulted in reduced level of genetic diversity in upstream populations. Gyrodactylus thymalli has a higher genetic diversity compared with G. salaris which has been introduced into new rivers (Hansen et al. Reference Hansen, Bakke and Bachmann2007b ; Anttila et al. Reference Anttila, Romakkaniemi, Kuusela and Koski2008). Low genetic variation may be an indication of newly introduced Gyrodactylus species. The results from this study should be taken into account when managing fish hosts and their parasites. Our study also indicates that sample size needs to be increased in such studies.
SUPPLEMENTARY MATERIAL
To view supplementary material for this article, please visit http://dx.doi.org/10.1017/S003118201500133X
ACKNOWLEDGEMENTS
We are grateful to Ole Haakon Heier, Jan Teigen, Erik Lien, Finn Gregersen, Henning Pavels and Camilla Fossum Pettersen for their assistance during sampling of grayling. Nanna Winger Steen and Emelita Rivera Nerli are acknowledged for their technical assistance in the laboratory and Kjartan Østbye for helpful discussions.
FINANCIAL SUPPORT
This research received no specific grant from any funding agency, commercial or not-for-profit sectors.
